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Introduce Ion 


The  multivariate  analyeia  of  heights  and  winds  is  one  of  the  key 
ingredients  of  any  global  data  assimilation  system  (GDAS),  e.g.,  McPherson  et 
al<  (1979),  Lorenc  (1981).  The  AFGL  GDAS  has  been  under  development  at 
APGL/LYP  for  several  years*  In  many  respects,  it  parallels  the  NMC  GDAS.  A 
GDAS  contains  three  main  components:  a  forecast  model,  an  analysis  method  and 
an  initialization  procedure.  discussefonly  one  part  of  the  analysis 

method,  namely,  the  multivariate  statistical  analysis  of  height  and  wind. 
Currently,  the  AFGL  GDAS  analyzes  moisture  in  a  separate  program  and  does  not 
analyze  surface  or  sea  level  pressure. 

This  documentation  does  not  attempt  to  provide  an  overview  of  the  entire 


AFGL  GDAS  or  even  of  the  entire  height-wind  analysis.  These  are  provided  by 
Norquist  (1986b).  -Ourcpurpose  jhere  is  to  describe  the  detailed  workings  of 


the  height-wind  analysis  program  ASAP1.  ft 
The  notation  used  in  this  documentat 
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gram.  Identification  between  the  , 


The  notation  used  in  this  documentation  does  not  directly  mirroT  the  <  , 
FORTRAN  symbolic  names  used  in  the  program.  Identification  between  the 
documentation  notation  and  the  symbolic  names  of  global  variables  (i.e.,  those*^ 
in  common  blocks)  is  made  in  Section  2.  Other  key  symbolic  names  are 
identified  or  defined  as  they  are  used.  Throughout  sections  2  through  A,  I 
use  the  following  notation: 

The  letters  u,  v  and  z  will  denote,  respectively,  the  eastward  wind 
component  (m/s),  the  northward  wind  component  (m/s)  and  the  height  of  a 
pressure  or  sigma  surface  (m).  Upper-case  letters  (U,  V,  Z)  without  any 
circumflex  are  reserved  for  observations.  Upper-case  letter  with  the  (~) 
circumflex  (U,  V,  Z)  are  reserved  for  analysis  or  first  guess  (i.e*,  fore¬ 
cast)  values.  The  residuals  are  identified  by  lower-case  letters  (u,  v,  z). 

For  example,  u  “  U-U.  The  letters  l  and  k  will  be  reserved  for  vertical 
indices.  The  letter  k  is  preferred  for  sigma  layer  variables,  e.g.,  is  the 
observed- temperature  (presumably  interpolated)  at  sigma  layer  k.  To 
distinguish  quantities  at  sigma  layers  from  those  at  sigma  levels,  I  will  put 
a  circumflex  (*)  over  the  first  character  of  the  subscript  of  a  sigma  level 
quantity,  e.g.,  Tj  is  the  first  guess  temperature  at  sigma  level  k.  There  are 
two  reasonable  ways  of  numbering  the  sigma  levels;  unfortunately,  both  are 
used  in  ASAP1.  First  of  all,  all  vertical  Indices  increase  with  increasing 
altitude  (and  decreasing  pressure).  Usually,  the  pressure  at  the  kch  sigma 
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level  is  Che  pressure  at  Che  bottom  of  Che  kth  signs  layer.  Thus,  F«  -  P* 
the  forecast  surface  pressure.  For  height  Z,  che  indices  are  usually  offset 
by  one;  che  height  ac  sigma  level  k  is  Che  height  of  the  Cop  of  Che  kch  sigma 

rw 

layer.  Thus  Zq  »  Z*,  the  copographic  height.  Sigma  level  usage  will  be 
defined  as  it  is  used  throughout  this  documentation.  The  letter  £  will  be 
preferred  for  pressure  levels  or  layers.  Usually,  data  is  associated  with 
pressure  levels;  usually,  a  circumflex  will  not  be  used  to  distinguish 
pressure  levels  from  layers.  Lower-case  letters  used  as  Indices  will  normally 
range  from  1  to  some  maximum  value,  which  will  be  denoted  by  the  corresponding 
upper-case  letter.  Thus,  Tj.*  k-l,K  is  the  set  of  first  guess  temperatures  at 
the  sigma  layers.  Pressures  associated  with  observations  will  be  denoted  P, 
pressures  associated  with  the  first  gusss  of  the  observations  evaluated  in 
sigma  coordinates  will  be  denoted  P,  pressures  associated  with  a  residual  will 

fW 

be  denoted  p,  and  mandatory  levels  will  be  denoted  P, .  For  example,  first 
guess  layer  pressures  are  P^. 

In  this  documentation,  an  observation  is  the  collection  of  all  data 
observed  at  a  particular  (X,$)  location.  Data  or  residuals  are  individual 
pieces  of  information.  Except  during  the  actual  update,  a  (u,v)  horizontal 
wind  is  considered  a  single  piece  of  information. 

Tables  1  and  2  summarize  the  above  discussion  and  list  other  commonly 
used  notation. 
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Table  1.  Variables  and  their  definitions 


Variable  Definition 


<U,V) 

z 

p 

~  a* 

(U,V),  Z,  P 


(u,v),z 

P 

T 

R 

Q 

iv  rv  rsj 

T,  R,  Q 
t,  r,  q 
X 
X 

$ 

a 

E(X) 

EP 

Eo 

B 

dsi 

£*££. 

Q(X)  * 

Q^X) 

Ql<x) 


p 


observed  horizontal  wind  (m/a) 

observed  height  of  pressure  or  signs  surface  (a) 

observed  pressure  (i*e*,  pressure  of  an  observation  before 
interpolation)  (ab) 

forecast  (first  guess)  or  analysed  values*  P  asy  also  denote 
the  pressure  of  an  Interpolated  observation.  These 
variables  refer  to  observation  locations  while  the  data  is 
read  and  prepared,  but  to  analysis  locations  during  the 
actual  update* 

residuals  of  (U,V),Z 

pressure  associate  with  residuals 

temperature  (Kelvin  or  Celsius) 

relative  humidity  (percent) 

specific  numidlty 

first  guess  values  of  (T,  R,  Q) 

residuals  of  (T,  R,  Q) 

generic  variable,  which  may  be  any  one  of  the  above 

longitude 

latitude 

sigma,  the  normalised  vertical  coordinate  defined  by  p/p* 

standard  deviation  of  X 

prediction  error 

observation  error 

buddy  check  box 

data  source  index 

data  type  (dsi/ 10) 

quality  mark  associated  with  X 

quality  indicator  associated  with  X 

quality  level  associated  with  X 

Qj  are  the  quality  control  flags  in  the  Level  II  data  sets.  Q 
are  the  generic  quality  marks  translated  from  the  Qj  which  are 
used  in  data  selection.  are  the  adjusted  Q  which  are  used 
in  the  buddy  check  flaggi.ig  procedure.  (See  Section  5.2*) 

horizontal  correlation 


v  vertical  correlation 


p  total  correlation 

e  normalised  observation  error 

P  mandatory  pressure  (mb) 


Table  2.  Subscripts  and  their  usage 


Subscripts 

Usage 

i,  J,  n,  n 

generic,  usually  used  for  observation  indices  and  pointers 

k,  i 

layer/level  indices*  k  la  preferred  for  o  layers /leva la, 
t  is  preferred  for  p  layers/levels*  The  (*)  circumflex 
denotes  a  sigma  level* 

n 

the  underground  layer/level 

* 

the  surface 

g 

the  particular  grid  point  under  consideration. 

(Either  (X,$)  or  (X,+,k).)  • 

V  £b 

layer/level  above,  below 

lu*  £t 

upper,  lower  layer/level 
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1 •  Global  Variables 


Tha  variables  scored  In  coanon  blocks  are  described  here,  cowson  block  by 
»oason  block,  in  alphabetical  order*  If  aore  than  one  syabollc  nsec  is  listed 
on  one  line,  then  the  different  nases  are  used  In  different  subroutines. 


Names 


AS  140  (common  block) 


Contends  Location*  and  asaodatad  variables  of  observations  and  tha 

currant  grid  point,  and  tha  Hat  of  observation*  selected  for 
the  analysis  at  tha  current  grid  point 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  A2R,  1986 

Referenced  by:  ASAP1,  ASAP 2 


Initialised  by: 


Symbolic  name 


Mean! 


YLG  4^,  grid  point  latitude  (degrees) 

XLG  Xg,  grid  point  longitude  (degrees) 

MF  mg,  map  factor 

XGG  Xg,  grid  point  coordinate,  stereographic  projection 

YGG  yg,  grid  point  coordinate,  stereographic  projection 

FIHG  +  (or  -)  1  for  a  grid  point  In  the  Northern  (or  Southern) 

Hemisphere 

NGBX,NGB0X  Ng,  number  of  observations  in  a  box  or  selected  for  analysis 

NAKEY(i)  list  of  observations  selected,  1*1,  Ng,  <  NSNDS 

NWORDS  1  +  5K  +  21^,  second  dimension  in  OBD  array  in  /RESID/ 

QLAT(n)  list  of  4n,  observation  latitudes 

QLON(n)  list  of  XQ,  observation  longitudes 

JDSI(n)  list  of  dsl0,  observation  data  source  Indices 

ABYLG  |4gJ ,  absolute  value  of  grid  point  latitude  (degrees) 

YGRAD  4g,  grid  point  latitude  (radians) 

XGRAD  Xg,  grid  point  longitude  (radians) 


AS  IAS  2  (coi 


!>lock) 


hirpoat!  Contains  constants  and  variables  for  ths  nod* Is  of  tho 

statistical  quantities*  Other  constants  are  hardwired  in 
SSOBBR. 


Author: 

Docuoentanlonx 

Kef are need  by: 

Initialised  by: 

D.  Norquist,  SASC,  1980  -  1986 

R*  Hoffman,  AER,  1986 

ASAP1 ,  ASAP2,  FLAGS 

AS  API 

Symbolic  nans 

Meaning 

PH(i) 

P^  ,  mandatory  pressures,  t*l,  12  (1000,  850, *..50  mb) 

ZBR 

limit  on  s  corrections  (250  m) 

UER 

limt  on  wind  corrections  (25  m/s) 

KH(i) 

k^»  horizontal  forecast  error  correlation  constants, 
i  -  1,3  ((2,  1.5,  1)  x  E-12m"2) 

KP 

k_,  vertical  forecast  correlation  constant  (5*0) 

DRAD 

»/4 

ERAD 

7e/4 

KHO 

kho»  horizontal  observational  error  correlation  constant 
(11.3  E-12m”2) 

RPO 

kpQ,  vertical  observational  error  correlation  constant  (8.3) 

KPOW 

k  -  ,  vertical  observational  error  correlation  constant  for 

Pv 

winds  (.’»*2) 

OMEGA 

Q,  rotation  rate  of  the  earth  (7.29116E-5  sec”*) 

NUKPOS 

look  up  table  for  satellite  height  vertical  observational 

error  correlations  (Norquist  1986b,  Table  3) 

EPZ 

Ep(ls),  EPE  (estimated  prediction  error)  for  heights  at 
mandatory  layers  for  this  latitude  (also  used  for  EaQ(Z^)) 

EPU 

E  (TJ  ) »  EPE  for  u  wind  component  at  mandatory  layers  for 
p  t 

this  latitude 
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E  (V.),  EPS  for  v  wind  component  at  mandatory  layers  for 
P  * 

this  latitude 

•  1  if  thia  la  the  flrat  (or  0  If  thla  la  not  the  flrat) 
analysis  in  a  sequence 

e1T^)»  P^^ictAon  error  growth  for  a  6  h  period 
at  mandatory  levels 

E„(^)t  analysis  error  limits  at  mandatory  levela 
Ef(^),  tropical  wind  EPE  at  mandatory  levela 


CCONST  (eoMon  block) 


VQtpOMS 

Author: 

Documentation: 
kef ere need  by: 
Initialised  by: 


Contains  constants  describing  the  analysis  and  fine  mesh 
grids  as  well  as  some  other  dimensional  constants* 

D.  Norquist,  SASC,  1980  -  1986 

R.  Hoffman,  AER,  1986 

AS API,  ASAP2,  PG,  FLAGS,  SETFG 

ASAP1 


Symbolic  name  Meaning 


MG 

NG 

HXG 

NYG 

A 

MAHDLVL 

NSNDS 

NLRHS 

DLAT 

DLON 


number  of  longitudes  in  fine  mesh  grid  (360) 
number  of  latitudes  in  fine  mesh  grid  (181) 
number  of  longitudes  in  analysis  grid  (61) 
number  of  latitudes  in  analysis  grid  (62) 

Ae,  radius  of  earth  (6.371  E6  m) 
number  of  mandatory  pressure  levels  (12) 

Maax,  maximum  number  of  observations  to  select  (8) 
maximum  number  of  pieces  of  data  to  select  (10) 

A+,  latitude  increment  for  the  fine  mesh  grid  (1°) 
AX,  longitude  increment  for  the  fine  mesh  grid  (1°) 
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Name:  DCONSf  (coaaon  block) 

Purpose:  Contains  physical  constants ,  dimensional  constants  and  the 

slsslng  value  indicator. 

Author:  D.  Norqulst,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 

Referenced  by:  ASAP1,  ASAP2,  CALCRES,  CNTOBS,  FG,  FLAGS,  LOWTMP,  LOWTMPS, 

MASTOR1,  MASTOR2,  MAST0H4,  MASTOR6,  POINTS,  PTOSIG,  SATLTHP, 
SETFG 


Initialized  by:  ASAP1 


Symbolic  name 


Mean! 


ML.IML 

PI 

DNN 

G,GRAV 

R,RG 

CP 

RKAPPA 

RK1 

MLP1.KPP1 

EPS 

MLRH 

LLL 


K,  number  of  model/analysis  sigma  layers  (12) 

«  (3.1415926536) 

missing  value  indicator  (-999.9) 

g,  acceleration  of  gravity  (9.8  m/s) 

R,  gas  constant  for  air  (287.05  J  kg"*  K"*) 

Cpt  specific  heat  of  air  (1005  J  kg"*  K"*) 

R/Cp 

R/Cp  +  1 

K  +  1,  number  of  sigma  levels  (13) 

molecular  weight  of  alr/molecular  weight  of  water  (1/0.622) 
Kr,  number  nf  sigma  layers  at  which  humidity  is  analysed 
(7) 

1  +  5K,  number  of  words  in  a  dry  profile  (see  /RESID/) 
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DEL  (coMon  block) 


Ruti 

Purpose:  Contain*  tha  residuals  for  single  level  data  as  calculated  by 

:  '  .  CALCRBS. 

Author I  D.  Norquist,  SASC,  1980  -  1986 

Bocuaantation:  R.  Hoff nan,  AE&,  1986 

Referenced  by:  CALCPJRS,  MASTOR2,  MAST0R6  . 

Initialised  by:  none 


Synbollc  naee 

Meaning 

DELT 

t. 

tenperature  residual 

DELU 

«» 

u*Trind  coeponent  residual 

DSLV 

v. 

v-wind  coaponent  residual 

DBLZ 

*» 

height  residual 

Name: 


FGDATA  (common  block) 


Purpose:  Contains  valuta  of  the  first  guess  fields  evaluated  at  the 

current  observation  location*  These  values  are  obtained  by 
interpolating  a  fine  aesh  (1*  x  1*)  grid  blllnearly  in  1  and 
4  (see  subroutine  FG). 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 

Referenced  by:  CALCRES,  FG,  MAST0R1 ,  MAST0R2,  MAST0R4,  HAST0R6 ,  PTOSIG, 

SATLTMP 

Initialized  by:  none 


Symbolic  name 


Mean! 


Z0(k)  Z£,  heights  at  oj^ 

eve 

T0(k)  Tk,  temperatures  at 

eve 

U0(k)  U^,  u  wind  component  at 

V0(k)  Vk,  v  wind  component  at  ok 

RHO(k)  R^,  relative  humidity  at  o^ 


FCFLDS  (common  block) 


Purpose  : 

Contains  two  latitudes  of  the  I*  x  1*  fine  mesh  grid  used  to 
evaluate  the  first  guess  at  the  observation  locations. 

Author: 

D.  Norquist,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman,  AER,  1986 

Referenced  by: 

FG,  SETFG 

initialised  by: 

none 

Symbolic  name 

Meaning 

TG  <i,J,k) 

UC  <i,j,k) 

VG  <i,J,k) 
RHG  <l,J,k) 
PSG  <l,j) 

ZSG  (i,j) 

** 

T(i, j ,k),  temperature,  k-l,K 

U(i,J,k),  u-wlnd  component,  k«l,K 

V(i , j ,k),  wind  component,  k-l,K 

R(i,j,k),  relative  humidity,  k-l,K 

a#  r 

P*(1,J),  surface  pressure 

ZA(i,j),  model  topography 

all  of  the  above,  lal,««.,360  Indexes  the  longitude  and 

j*l,2  Indexes  the  two  latitudes  currently  in  memory. 
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Name: 

1FLG  (common  block) 

Purpose: 

Contains  flag  arrays  used  in  subroutine  FLAGS  to  perform 
buddy  check.  Maximum  number  of  observations  to  be  checked  in 
a  box  is  currently  100. 

Author: 

D.  Norquist,  SASC,  1980  —  1986 

Documentation: 

R.  Hoffsan,  AER,  1986 

Referenced  by: 

FLAGS 

Initialized  by: 

none 

Contents: 

Symbolic  name 

Meaning 

IFZTI 

*Fj 

IFZTJ 

,p! 

IFWI 

“Fi 

IFWJ 

Vpi 

3 

See  FLAGS(O)  for  a  detailed  discussion. 


1A 


Kane: 

LOUT  (common  blocV.) 

Purpose: 

Contains  constants  nesded  to  apply  Flattery  algorithm  to 
convert  layer-layer  temperatures  to  sigma  layer  temperatures. 

Author: 

D*  Norquiet,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman,  AER,  1986 

Referenced  by: 

LOWTMP 

Initialised  by: 

ASAP1 

Contents : 

• 

Symbolic  name 

Meaning 

AL,  AC 

WL 

BL,  BC 

«U 

AT 

at 

ATA 

(A1*)’1 

Refer  to  Appendix  A  of  Norquiet  (1986b,  pp.  55-58)  for  details* 
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Name : 

RES1D  (common  block) 

Purpose: 

Contain*  the  residuals  calculated  by  the  MASTOtn  routines* 

Author: 

D.  Norquist,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman,  ASR,  1986 

Referenced  by: 

ASAP1,  ASAP2,  CONTOBS,  FLAGS,  HAS TORI,  MASTOR2,  MASTOR4 , 
MASTOR6 ,  POINTS 

Initialized  by: 

none 

Contents: 

Symbolic  name 

Meaning 

OBD  (n,  1) 

contains  Information  for  obaervatlon  n- 1,6000.  The  various 

residuals  are  Indexed  by  t.  They  are  stored  in  the 

order: 

P*.  l*£,  Q<*£)»  \»  vk»  ^uk»vic>»  k-1,8), 
lrk»  Q(rk),  k-l,KrJ.  Thus  8-1,  1+58+21^  -  1,75. 
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SICK  (coMon  block) 


Mom: 

Purpose:  Contain*  variable*  describing  the  vortical  structure  of  the 

model* 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  ABR,  1986 

Referenced  by:  ASAP1,  ASAP2,  FG,  FLAGS 
Initialised  by:  no»u 

Symbolic  name _ Meaning _ 


SIGMA(k),  SL(k) 

SLU 

SI(k) 


o^,  k«l,K,  sigma  in  layer  k 
au,  sigma  in  underground  layer 

o£,  k»l,K+l,  sigma  at  Interface  at  bottom  of  layer  k 


Name:  UADATA  (common  block) 

Purpose:  Contains  data  as  road  In  by  MASTORn  routines  from  WE  Level 

II  data  sets.  Usage  varies  with  context.  E.g.f  ia  MAST0R1 
Zt  is  the  height,  at  Pftt  while  In  MA8XOR4  it  la  the  height 
relative  to  Zj  of  P^. 

Author:  D.  Norqulst,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 

Referenced  by:  MAST0R1 ,  HASTOR2,  MAST0R4 ,  MASTOR6,  PTOSIG 
Initialised  by:  none 


Symbolic  name 


Meaning 


PU> 

V 

z(t) 

V 

T(t) 

UU) 

V 

V(i) 

V 

Qd) 

V 

IQZ(i) 

Q(2j 

IQT(t) 

Q(T, 

IQW(t) 

Q(Uj 

<ui 

iQQU).IQCQ(i) 

Q(Q, 
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Mane: 

UASIGNA  (common  block) 

Purpose: 

Contain*  tha  observations  and  thalr  associated  quality  marks 
Intarpolatad  to  nodal  algua  layara/lavals. 

Author} 

D*  Norqulst,  SASC,  1980  -  1986 

Documentation: 

8.  Hoffman,  AER,  1986 

ioforoneod  by: 

CALC RES,  MAS TORI ,  MAST0R2 ,  HAST0R4 ,  MAST0R6 ,  PTOSIG 

Initialised  by: 

P 

1 

Symbolic  name 

Meaning 

ZBO(k) 

Z*,  k«l,K,  height  at  o*+1 

TBO(k) 

Tk,  k-1 ,K,  teaperatura  at 

UBO(k) 

Ukt  k-1 ,K,  u -component  mind  at 

VBO(k) 

k-1  ,K,  v-componant  wind  at  o^ 

QBO(k) 

Q^,  k-1  ,Rr,  specific  tumidity  at  0* 

IQZO(k) 

Q(Z* ),  k-l,K,  quality  nark  asaodatad  with  z* 

IQTO(k) 

Q(Tk),  k-1  ,K,  quality  mark,  associated  with  Tk 

IQWO(k) 

Q(,Jk,vk.)»  h-l.K,  quality  mark  associated  with  (U^,^) 

IQQ0(k),IQCQ0(k) 

0(0^),  k-ljKy,  quality  nark  associated  with 
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3.  Progran  Units 


Ths  asln  progran  AS API  and  tha  various  aubprograaa  It  Involvaa  ara 
daaerlbad  bars  In  alphabatlcal  ordar>  Within  aach  description,  thara  ara 
various  nuabarad  subtactions.  Tha  start  o 2  aach  of  thasa  blocks  Is  Identified 
tar  a  coaaant  card  In  tha  FORTRAN  coda. 


ASAP1  (sain  program) 


Name : 

Purpose: 


Author: 

Documentation: 

Referenced  by: 

References: 

Commons  used: 

Arguments: 

I/O  units: 
Description: 


Performs  the  statistical  analysis  procedure  following  tergman 
(1979)  and  Day  and  Morone  (1985)  as  described  in  Norquiat 
(1982,  1983,  1984,  1986a,  1986b)  and  Halberatam  at  al« 

(1984).  The  moisture  analysis  is  now  separate  but  moisture 
residuals  are  calculated  here  for  the  GWE  Level  II  data  and 
stored  with  their  associated  Identification  and  location 
parameters  on  units  8  and  3.  AS API  calls  several  routines  to 
read  in  the  data,  interpolate  to  sigma  coordinates,  and 
perform  quality  control.  It  then  loops  over  all  the  analysis 
grid  points,  oelecting  observations  and  calling  A8AP2  to 
actually  perform  the  analysis.  Mote  then  an  observation  may 
be  an  entire  RAOB  report  or  a  single  level  AIREP.  The  KA0B 
contains  many  pieces  of  data,  the  AIREP  only  one. 

D.  Ncrquist,  SASC,  1980  -  1986 

R.  Hoffman,  AER,  1986 

None 

ASAP2,  CNTOBS,  FLAGS,  MAST0R1 ,  MAST0R2 ,  MAST0R4,  MASTC16, 
POINTS,  SETFG. 

AS140,  AS1AS2,  CCONST,  DCONST,  RESID,  SIGK 
None 

2,  3,  7,  8,  10,  INPUT,  OUTPUT 


(1)  Initialize  constants 


AS  API 


(2)  Bead  In  na Delist  containing  variables  describing  date  and  tine  and 
whether  this  is  the  first  analysis  in  a  sequence. 

(3)  Open  unit  2,  the  unpacked  GWE  Level  £1  data  set,  described  by  Norqulst 
(1984,  p.  17*33)  and  check  header  record  date/tlae  against  that  read  via 
naaelist  (2). 

(4)  Set  up  the  vertical  structure  in  /SICK/.  (See  Section  5.4.) 

(5)  Read  in  type  1  data  (RA0B5)  from  file  7  on  unit  2  (MASTOR1).  In  (5)  * 

(8)  the  call  to  SBTFG  initialises  the  horisontal  interpolation  procedure. 

(6)  Bead  in  type  4  data  (SAIEMS)  fron  file  5  on  unit  2  (MAST0R4). 

(7)  Read  in  type  2  data  (A1REPS)  fron  file  3  on  unit  2  (MAST0R2). 

(8)  Read  in  type  6a  data  (SATWINDS)  fron  file  7  on  unit  2  (MAST0R6). 

(9)  Count  valid  data  read  in  and  display  totals  (CNTOBS). 

(10)  Gross  and  bu  'dy  check  data  (FLAGS),  assign  points  to  each  observation 
(POINTS).  These  points  will  be  used  in  the  observation  selection 
procedure.  Count  renalnlng  data  and  display  totals  (CNTOBS). 

(11)  Store  residuals  on  unit  8  and  identification/location  parameters  on 
unit  3. 

(12)  Beginning  of  analysis  phase  ((12)  *  (27)].  Select  tropospheric  value  of 
kh  (2.08*12  n~^)  given  by  Dey  and  Morone  (1985).  (See  (23).)  Begin  loop 
over  all  latitudes  in  the  analysis  grid. 

(13)  Deternine  Che  latitude  of  the  current  grid  row.  The  sine  of  $g  is 
read  fron  unit  7  if  Guassian  latitudes  are  used  (ITLAT  -  1),  otherwise 
regularly  spaced  $g  are  used.  (Currently  ITLAT  io  hard-wired  to  1.) 


ASAP1 


(14)  If  this  is  the  first  analysis  in  a  sequence  read  in  the  standard 

estimated  prediction  errors  (EPE)  at  the  mandatory  levels  from  unit  10* 

In  this  case  the  EPE,  E  (2.),  E  (tJ. ),  E  )  do  not  depend  on  longitude* 

pi  PI  Pi 

(15)  Begin  loop  over  all  longitudes*  Define  Xg* 

(16)  If  }$g|  >  70,  convert  (Xg,  ♦g)  to  (xg,  yg)  the  coordinates  in  a  polar 
stereographic  projection  true  at  the  pole*  (See  Section  5*5*) 

(17)  Define  the  latitude /longitude  observation  selection  box*  The  box  is 
±  15°  in  latitude  and  ±  30*  in  longitude  if  |^-|  <  60**  (Longitude 

— 9  ® 

window  would  be  a  AO*  if  kjj  <  1*8  E  -  12  m  »  but  this  trill  not  occur  for 
current  parameter  settings*)  For  j+gl  >  60*'  all  longitudes  trill  be 
considered.  (See  (20)  below.) 

(18)  Begin  data  selection  phase  1(18)  -  (25)].  (See  Norqulst  1986b,  pp.  9- 

10.)  Set  to  0.1*  Start  loop  on  data  type*  NTY  -  1,  2,  3,  A,  which 

min 

correspond  to  type  -  1,  A,  2,  6a*  If  (i.e.  NSNDS)  observations  have 

been  selected  exit  loop.  Begin  loop  on  observations  within  this  type* 
This  algorithm  implies  that  if  *»ax  type  1  observations  are  found  with 
uzz  >  u^n  other  data  type3  are  not  considered,  but  all  the  type  1 
observations  are  considered*  Similarly  for  type  1  and  A  combined,  type 
1,  4,  and  2,  combined,  etc. 

(19)  Check  that  $n  la  within  the  selection  box.  If  not  make  use  of  fact  that 
data  are  sorted  by  latitude.  If  1®  south  of  box  go  to  end  of 
observation  loop.  If  $n  is  north  of  box  go  to  end  of  type  loop. 

(20)  If  |$g|  <  60,  check  that  An  within  the  selection  box.  If  not  go  to 
end  of  observation  loop*  Note  the  check  when  the  box  straddles 
Greenwich,  i.e.  when  XMX  <  X^p.  In  this  case  XMX  <  Xn  <  X^u  implies 
Xn  is  not  in  the  box. 

(21)  If  |$gj  >  70  then  transform  observation  location  (Xn»  An)  to  polar 
stereographic  projection  coordinates  (xn,  yn)  and  calculate  dgn  the 
distance  between  the  grid  point  and  the  observation. 


2A 


A5AP1  i 


(22)  Else  If  |+g|  <  70  Chen  calculate  distance  with  approximate  great  circle 
foraula  (Schlatter,  1975). 

dgn  “  *2[<C08  <♦>  ^>2  +  <A*>2] 

(23)  Calculate  horizontal  correlation  pZ*  »  exp  (-k.d  ^).  If  the  number  of 

gn  h  gn 

observations  selected,  H  <  HBax  skip  observations  with  p  <  umln*  II 
M  >  skip  observations  with  p*  <  P*Bin«  Here  *  Indicates  that  p  has 

been  multiplied  by  the  number  of  points  assigned  to  the  observation  in 

(10). 

(24)  Insert  the  current  observation  in  the  list  of  observations  for  this 

analysis  point  stored  in  descending  p*  sort  order.  If  M  was  already 
Mmx,  drop  the  last  one  in  the  list  and  reset  y  If  M  was  Mmax  -  1, 

set  p*Bin  to  the  last  p*  in  the  list.  The  list  is  maintained  in  two 
arrays:  NAKEY(m)  -  n,  which  points  to  the  observation  arrays  and  MUTTXY 
(a)  -  P*n« 

(25)  End  of  data  selection  loops  on  observation  and  type  (18). 

(26)  Apply  statistical  Interpolation  procedure.  (ASAP2). 

(27)  End  of  loops  over  analysis  longitudes  (15)  and  latitudes  (12). 

(28)  End  of  main  module. 
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ASAP2 

Name : 

Purpose: 


Author: 

Documentation: 

Referenced  by: 

References : 

Commons  used: 

Arguments: 


ASAP2  (subroutine) 

Performs  update  of  all  variables  at  all  sigma  layers 

for  a  single  8r*d  location.  The  procedure  follows 

that  of  Dey  and  Morone  (1985)  and  Bergman  (1979)  as  described 
by  Norquist  (1986b,  pp.  10-18). 

In  brief,  ASAP2  obtains  (a)  the  first  guess  values  and 

(b)  the  expected  prediction  errors  at  the  analysis  ppint. 
Then,  for  each  layer,  including  a  subsurface  layer,  ASAP2 

(c)  selects  observations  highly  correlated  with  the  grid 
point  variables,  (d)  calculates  the  observation  location-grid 
point  prediction  error  correlations  for  the  selected 
observations  (i.e.,  the  r.h.s.  of  the  normal  equation), 

(e)  calculates  the  matrix  of  observation  location-observation 
location  prediction  and  observational  error  correlations  (the 
same  matrix  is  used  for  2,  B  .and  ^),  (f)  solves  the  normal 
equations  for  the  nondimensional  analysis  weights  and 
(g)  calculates  the  corrections  and  expected  analysis  errors. 
Finally,  (h)  z  corrections  are  converted  to  T  corrections, 
and  the  first  guess  is  updated  and  stored  along  with  the 
corrections  and  estimated  analysis  errors. 

D.  Norquist,  SASC,  1980  -  1986 

R.  Hoffman,  AER,  1986 

ASAP1 

CHLSKY,  ESCBER,  LOWTMP,  RECALC 

AS140,  AS1AS2,  CCONST,  DCONST,  RESID,  SIGK 

1GP  Ag  index 

JGP  dg  index 

IKON  month  of  the  analysis 


26 


ASAP  2 


I/O  unless  *»  6,  10,  11,  OUTPUT 

Description:  Subscript  snd  pointer  usage  In  ASAP2  are  sore  complicated 

than  elsewhere: 


Symbolic  name  Meaning 


IN,  JN 


1,  J  Index  piece  of  data  which  is  used  in  analysis 


NVI,  NVJ  0^(1)  data  type  (1,2,3)  -  (s,u,v)  or  (1,2) 


NN 


m(i)  intermediate  pointer  to,  n  and  1  once 


NI,  NJ 

n(i)  -  n(m(l),DT(i)) 

LI,  LJ 

£(i)  -  tOKD.Djd)) 

ZB(NI,LI) ,  etc. 

x(i)  -  x(t(i),n(l),D 

-  (a , (u , v) ) 
Dj(l)  is  known 
observation  index  within  local  array 
layer/level  index 

the  1th  datum* 
and  v 


It#  is  stored  in  one  of  three  arrays  for  z,  u 
separately,  depending  on  the  value  of  DT(i) 


T(i))  is  the  value  of 


PH(NI,LI),  etc.  p(l)  -  p(l(i),n(l),DT(i))  is  the  value  of  pressure 

associated  with  the  1th  datum*  It  is  stored  in  one  of  two 
arrays  for  s  and  (u,v)  separately,  depending  on  the  value  of 
Dr(i) 


E 


expected  standard  deviation,  e*g. , 
£^(2^)  EPE  of  the  first  guess 

Ea(zkg)  EAE  (becomes  Ea0(2k)) 


EAE  of  previous  forecast 


ASAP2 


(a)  -  (1)  -  (2)  Obtain  the  first  guess  values 

(1)  From  unit  ■  4,  read  in  first  guess  values  ($k»  ^k»  ^*)  for  current 

location*  Compute  ?k  and  ?u*  Here  and  below  the  subscript  u  denotes  the 
underground  layer* 

(2)  Calculated  G  ,  the  coefficient  of  geostrophy  at  the  grid  point* 

6 

(See  Norquist  1986b,  p.  14  and  Section  5*6*) 

(b)  ”  (3)  -  (7)  Determine  the  expected  prediction  errors  (EPS)* 

(3)  Begin  determination  of  EPE  (expected  prediction  errors)*  If  this  is  the 
first  analysis  in  a  sequence  (IFRST  -  1),  then  use  the  a  priori  WIC 
values  Interpolated  linearly  in  in  p  to  the  sigma  layer  pressures  Pk* 
Above  50  mb,  the  50  mb  value -is  used,  below  1000  mb  extrapolation  is 
used.  Since  the  NMC  values  do  not  depend  on  Ag  they  are  read  outside  th? 
longitude  loop  in  ASAP1(14)  and  stored  in  /AS1AS2/*  The  EPE  values 
EpCZfc),  EpO^),  Ep(Vk)  are  stored  in  variable  EPE.  Note  that  EpC\)  is 
the  EPE  of  Z  at  ok  not  at  o*+1  (see  (7)). 

(4)  Otherwise,  if  this  is  not  the  first  analysis  in  a  sequence  (IFRST  "0), 
then: 

From  unit  -  10,  read  EaQ(Zk),  the  expected  analysis  error  of  the  previous 
analysis.  Interpolate  E._(Z.  )  the  EPE  growth  during  the  interval  between 
analyses  from  the  mandatory  layer  values  of  Dey  and  Morone  (1985,  Table 
3,  p.  309).  Also  interpolate  E  (Z. ),  the  climatological  maximum  analysis 

•1C  M 

height  errors  from  Table  2  of  Dey  and  Morone,  and  E^(Uk),  the  tropical 
wind  EPE  from  Table  4  of  Dey  and  Morone.  In  all  three  cases  the 
interpolation  procedure  is  identical  to  that  used  in  (3).  Next,  force 
Efl0  to  be  in  the  interval  (E^E^).  Then  set  Ep(Zk)  •  Ea0(Zk)  +  EAT(Zk), 
and  calculate  Ep(Uk)  ■  Ep(V^)  using  equations  (2.12)  and  (B.ll)  of  Dey 
and  Morone  (1985)  as  necessary 

Ep(Uk)  -  Ep(Vk)  -  (g/|f|)(2kh)1/2  Ep(Zk) 

W  ■  VV  w  +  (1 


28 


ASAP2 


vhtra  la  sac  as  la  PLAGS(7).  The  blsndsd  vslass  Eg  ars  used  only  If 
|tg|  <  23*  (Hots  Chat  typographical  errors  In  (B. 11)  have  bean  corrected 
here*) 

(3)  Initialise  variables*  setting  all  corrections  (i^,  u^,  v^)  to  sero  and 
all  estlaated  analysla  errors  (EAE)  to  the  EPE,  e.g. ,  E#(R^g)  »  Ep(Zk). 
(6)  Loop  over  observations  previously  selected  in  ASAP1  (24)*  i.e.,  consider 
observations  in  the  list  HAKBT(l)*  1*1*NGB0X.  Extract  data  (pk*  u^,  v^, 
QCu^*  v^),  p£,  a£  and  q(Z£))  f roe  /RESID/  for  each  of  these  observations 
for  each  k  for  which  data  is  present  and  store  in  local  arrays.  Observa¬ 
tions  with  no  data  present  are  skipped.  Since  layers/levels  with  no  data 
are  skipped*  the  pressures  associated  with  the  data  are  stored.  That  is, 
uin  is  the  u  residual  at  the  tth  layer  which  had  wind  or  height  data 
present  in  the  nth  observation*  and  the  pressure  of  u£n  is  pln.  the 
nuaber  of  layer /levels  extracted  is  stored  in  LQ.  Local  arrays  of  dsin, 

a  ,  X  are  also  foraed.  If  |ip  |  >  70*  (u*v)  are  converted  to  polar 
n  n  g 

stereographic  corrdinates.  If  N*  the  nuaber  of  observations  extracted* 
is  sero*  go  to  (49). 

’(7)  Assign  the  EPE  for  the  residuals*  Sp(s£n)*  Ep(uln)  and  2p(vAa)  by  finding 
the  grid  point  signs  layer  pressure  P^  closest  to  p*Q  for  z  or  p^n  for 
(u*v)  and  assigning  the  corresponding  grid  point  EPE  to  the  observation 

*» .  *,(«*„>  -  v\,)- 

(8)  Initialize  IWFLG.  Begin  aain  loop  (8-46)  over  layers*  k»l*K+l.  Layer 

1  is  the  underground  layer.  Select  kh  as  in  FLAGS(7).  P  is  the 

*  <v  © 

insure  at  the  grid  point  (l.e.*  Pk,  or  Pu  if  k-K+1). 

(c)  ■  (9)-(14)  Selects  the  data  to  be  used  at  this  layer. 

(9)  Loni  over  all  observations  n  ■  1*N.  Set  L  ■  Ln« 

(10)  Calculate  d  *  the  distance  froa  the  observation  to  the  grid  point. 

F<*rUgl  >  70**  the  polar  stereographic  projection  is  used.  Calculate 

G_*  the  coefficient  of  geostrophy  at  the  observation,  and  y zx,  the  zz 

8 

horizontal  correlation.  (See  Section  5.6  and  Dey  and  Morone  (1985, 

Eq.  2.8) . ) 


I 
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(11)  Begin  loop  over  all  layer/levels  In  Che  observation  £  -  1,L.  If  a  a 

residual  is  present,  then  add  n  and  £  to  list  of  potential  s  observation 
and  levels  (variables  NZT(a),  LZT(a),  a-l,NTOTZ).  Calculate  the  vertical 
correlation  between  p£n  and  Pg  using  Eq.  (2.9)  of  Dey  and  Morone 
(1985).  Obtain  the  estimated  observation  error  (EOE)  for  this  residual, 
V'in*  <ES0BER>*  Define  e>  a8  E0(aJn)/Ep(a£n).  Calculate  p^*,  p^u 
and  P^|v»  the  adjusted  total  correlations  between  the  height  residual  and 
the  grid  point  variables.  (Note  that  a  implies  a  particular  £  ard  n, 
while  g  implies  k.)  These  adjusted  correlations  are  the  correlations 
obtained  following  Dey  and  Morone  (1985,  Eqa.  (2.7),  (A2),  (A3),  (B2), 

(B3)  as  appropriate),  divided  by  (1  +  e  *)•  These  adjusted  correlations 

a 

are  the  normalised  weights  one  would  obtain  for  updates  made  using  each 
residual  separately.  The  auta  of  the  absolute  values  of  these  weights  is 


P 


Ip,,ui 

,Hag 


+  h 


•*v 

■8 


I 


For  the  underground  layer  (i.e.,  if  k-K+1),  p(D  ■  |p'*x|  • 

(12)  Repeat  (11)  for  wind  data.  If  wind  data  are  present,  then: 

Add  n  and  £  to  list  of  potential  (u,v)  observations  and  layers  (variables 
NW(m),  LW(m),  m-l,NTOTW).  Calculate  vertical  correlation  between 
p£n  and  pg*  Obtain  E0^u£n^  *nd  define  em  as  Eo^u£n^Ep(u£n^*  ^Not*  that 
Ep<u£n)  is  used  for  v  as  well.)  Calculate  P^*»  Pi£*»  P^U»  PagV»  Pag1 
and  pmgV>  the  adJusted  correlations,  as  in  (11),  but  here  using  Eqs  (A4)- 
(A9)  or  (B4)-(B9)  of  Dey  and  Morone  (1985).  Set 


(2) 


m 


[I. 


,'UZ| 

mg  1 


vs  i 


i +  i»n  ♦  ip^'i 


UVi 


vu 


mg 


|pig 


I  + 


vv 


K 


l]/2 


For  the  underground  layer,  set  p£^  ■  [|p,ur| 
End  of  loop  (9). 


(13)  Determine  the  NLRHS  (currently  10)  largest  |p|  greater  than  0.1  from  the 

combined  set  of  and  p^.  At  this  point,  a  (u,v)  wind^is  considered 

D  D 

a  single  piece  of  data.  Pointers  to  the  selected  residuals  are  saved  in 
m(i),  DT(si),  i-l,Y  (variables  NNRHS(i),  NVRHS(i),  i-l,NKT).  The  insert 
sort  algorithm  (see  ASAP1(24))  puts  these  pointers  in  the  list  in 
descending  |pj  order. 

(14)  If  1-0  and  k  <  K,  go  to  (45).  Else  if  1-0  and  k-K+1,  go  to  (46). 

Else  I>0;  continue  with: 
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(d)  »  (15)  -  (20)  Calculate  Che  observatlon-gridpoint  prediction  error 
correlations  for  Che  selected  residuals.  At  this  point,  the  u  and  v  wind 
components  begin  to  be  treated  separately. 


05) 

06) 

07) 


08) 


09) 


(20) 


Set  j-0.  Loop  over  selected  data  i*l,I  •  Extract  DT(i)  and  n(i)  and 
l(i)  froe  i(i) 

Repeat  (10),  calculating  d  ,  G_  and 

M  “O  U 

If  IXj.(l)  ■  1,  then  the  current  residual  is  a  height  residual:  Increment 

J.  Add  it  to  the  list  of  pointers,  by  setting  m(J)  ■  l(i)  nnd  0T(j)  *  1 

(variables  NNKHSV,  NVRHSV).  Repeat  (11)  but  store  p,M,  p,EU,  and  p’8v 

ng  ng  ng 

as  the  three  r.h.s.  of  the  normal  equations  (variable  (RHSV). 

If  DfO)*!,  then  the  current  residual  is  a  wind  residual: 

Add  the  u  component  to  the  list  of  pointers  as  in  (17).  Calculate 

P,u*»  P,uu»  and  p,uv  as  in  02).  Store  these  adjusted  correlations  as 

mg  mg  mg 

the  three  r.h.s.  of  the  normal  equations. 

If  Dj(i)-2,  repeat  (18)  for  the  v  component.  Note:  Dj(j)  "  3»  This  is 
the  point  where  a  wind  residual  is  treated  as  two  pieces  of  data. 

Calculate  p,v*,  p,vu,  and  p,vv  as  in  (12)  but  store  as  the  three  r.h.s. 

mg  mg 

of  the  normal  equations. 

End  of  Loop  (15).  Set  I-J.  The  r.h.s.  and  the  pointers  DT(i),  m(i)  are 
complete. 


(e)  ■  (21)-(32)  Calculate  the  matrix  of  observation-observation  correlations. 


(21)  Loop  on  all  pieces  of  data  i-l,I.  Extract  D-j(i)*  n(i),  £(i),  p(i)  and 
Bp(x(i)).  Obtain  E0(x(i))  (ESOBER)  and  set  s*  -  E0(x(i))/Ep(x(i)). 

(22)  Calculate  variables  needed  for  the  distance  calculation  and  the 
correlation  calculation  which  depend  on  the  ith  piece  of  data.  Calculate 
Gi»  the  coefficient  of  geostrophy. 

(23)  Loop  on  all  pieces  of  data  j«i,I.  If  j-i,  go  to  (31).  Otherwise,  obtain 
Dx(j),  o(j ),  l(j),  p(j),  Ep(x(j )),  E0(x( j ))  and  ej  as  in  (21). 

(24)  Calculate  the  distance  d^j,  the  coefficient  of  geostrophy  Gj,  and  other 
variables  needed  for  the  correlation  calculation. 

(25)  Calculate  variables  for  the  vertical  correlation  calculation.  Calculate 


the  horisontal  correlation  y**. 
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(26)  Set  horizontal  observation  error  correlation  (Morqulst,  1986b,  pp.  14-16) 
as 


as  „ 

M0ij 


1 

exp(-kho<l2) 

0 


If  use  variable  and  sane  observation 
(D*(j)  -  DjU)  And  n(J)  -  n(i» 
for  type  1  or  2 

for  satellite  heights 
otherwise 


(27)  Calculate  vertical  error  correlations.  Let  x  -  tn*(p(i)/p(j)) 

Define 

r 

for  heights,  type  1  or  2 


v 


OiJ 


O  +  hpQx)' 


ci  ♦  k 


Interpolate  linearly 
in  x  in  look-up  table 


for  winds,  type  1  or  2 


for  satellite  heights 
(Norquist,  1906b,  Table  3) 

otherwise 


and 


vij  *  <l  +  kp*> 


-l 


(28)  For  U  |  >  70,  calculate  p  and  n  using  formulas  froa  appendix  A  of  Dey 

w 

and  Horone  (1985).  Each  of  the  9  cases  ss,  su,  sv,  us,  uu,  uv,  vs,  vu, 
vv  are  handled  separately  here.  The  statistical  aodel  for  p  is  also  used 
for  n,  only  the  constants  and  kp  ara  by  k^,  kpQ  and  kj^. 

(29)  For  | 4  |  <  70,  calculate  p  and  n  using  formulas  froa  appendix  B  of  Dey 

© 

and  Horone  (1985). 

(30)  Set  matrix  entries  R^  -  R^  -  p^  +  CjSjHjj  *  Go  to  (32) 

(31)  Set  matrix  entry  R^  ■  1  ♦ 

(32)  End  of  loop  (23)  on  j.  End  of  loop  (21)  on  1. 
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(f)  •  (33)  -  (39)  Solve  Che  nor  Ml  equation* 


(33)  Begin  loop  on  data  type  Dj(g)  ■  1,3*  For  underground  layer  k  •  K+l,  aklp 
caaea  for  D^(g)  -  2,3,  since  only  the  a  analysis  Is  needed* 

(34)  Move  which  was  stored  In  ( 17)— ( 19)  In  variable  RHSV  to  the  vector  of 
rhs  b  (variable  RHS)* 

(35)  If  only  one  piece  of  data  is  present,  solve  for  this  special  case 
a'j  •  bj/Rjj.  Go  to  (38). 

(3d)  Solve  normal  equations  Ra*  -  b  using  Cholesky  algorithm  (CHLSKY). 

(37)  If  algorithm  falls,  then:  Eliminate  one  residual  which  is  highly 
correlated  with  another  residual  (RECALC).  If  only  one  residual  is  left, 
solve  for  this  special  case  as  in  (35);  otherwise,  apply  Cholesky 
algorithm  again  (CHLSKY).  If  Cholesky  algorithm  falls,  go  to  (37). 

(38)  Calculate  norMllsed  squared  analysis  error  Ag  from  Bergman's  (1979) 

Eq.  2.11.  If  any  of  the  **<  •>  1.1  (and  if  I>1),  go  to  (37),  since  this 
is  considered  a  failure. 

(39)  If  the  norMllsed  squared  analysis  error  Ag  is  less  than  sero  or  greater 
than  one,  skip  the  analysis  for  this  variable  by  setting  the  correction 
C-0,  Ag-0  and  1*0  and  going  to  (42a). 

(40)  Convert  a^  to  *t  the  dimension* Used  weights  using  (Norqulst,  1986b, 


Bq.  U) 


•i  -  **i  Ef/Ep(x(i)) 


where  E^  ■  EPB(k,DT(g))  for  k  <  K  and  Eg  •  Ep(Z|)  for  the  underground 

layer.  If  |*  j  <  10,  then  s  observations  are  not  allowed  to  affect  wind 

analyses;  set  a^-0  in  this  case.  If  <  10,  then  wind  observations 

are  not  allowed  to  affect  height  analyses;  set  a^-0  in  this  case.  The 

reasoning  is  that  if  U  I  <  10,  then  G  -  0  and  \i  -  u  -  0,  while  if 

,Tg'  g  Hsu  Hsv 

1 9^ |  (  10,  then  C^  •  0  and  *  y  -  0.  (See  Dey  and  Morone,  1985, 
Eqs.  (B2)-(B5)). 


(41)  Recalculate  Ag  as  in  (38),  but  skip  data  with  aj-0.  If  Ag  <  0  (or 
Ag  >  1),  print  a  measage  and  sot  Ag  to  0  (or  1). 

(42)  Calculate  C,  the  correction.  C  *  Za<x(i) 

i  1 

(43)  Perform  gross  error  check  of  C  against  ZER  (250  m)  or  WER  (25  m/9) 
(Norqulst,  1986b,  p.  17).  If  C  passes  the  check,  go  to  (44). 

(43A)  Set  C-0,  Ag-1.  If  C  is  a  wind,  set  the  flag  (IWFLG(K))  so  that  both  u 
and  v  corrections  are  set  to  zero  in  (48). 
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(44)  If  k-K+1,  set  Zug-C;  otherwise,  lit  B^s,^),  E^u^)  or  Ea(vkg)  equal  to 

aE/2  tiBea  Ep^>»  Ep<“k*  or  Bp^k  *  and  *kg*  \g  or  vkg  eo  C  if 

DTCg)  -  1,  2  or  3.  (So*  Eq.  2.11  of  Bergaan,  1979.)  Bad  of  loop  (33)  oa 

DT(g).  Go  to  (46) 

(45)  Set  *kg  ■  ukg  “  v^g  ■  0  and  aasodated  BAB  to  the  corresponding  BPE. 
(Llke(5)  but  for  single  k  only.) 

(46)  End  of  naln  loop  (8)  on  layers 

(47)  Convert  *kg  to  tkg  (LOUTMP)  and  update  first  guess  \  \  +  tkg. 

(48)  Update  (\»\)  by  adding  (Ukg»vkg);  first  converting  analysis  winds  back 
to  (A, 4)  coordinates  if  |4^|  >  70. 

(49)  Output  ?R,  llk,  ^k  to  unit  6  and  B#(akg)i  B>(ukj),  B#(vkj)  along  with 

*kg»  ukg*nd  vkg  to  unit  11#  Not*  that  Ba(*kg)  hBH  be  used  for  the 
next  analysis  (see  (4)).  Return. 
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Mimi  BILINR  (subroutine) 

flrpMtt  IntirpolitM  a  single  field  bill  nearly* 

Author i  D.  Norqulst ,  SASC,  19S0  -  1986 

Documentation:  R.  Hoffman,  AIR,  1986 

leferencad  by:  FG 

Raferencea:  none 

Commons  used:  none 

Arguments:  A  f ^ , 'field  array 

(X,T)  (x,y)  location  in  grid  units 

NX  first  A  array  dimension 

NT  (not  used) 

ANS  f(x,y),  interpolated  value 

1/0  units:  none 

Description: 

Determine  the  coordinates  (l,j)  of  the  grid  point  In  the  lower  left  hand 
corner  of  the  grid  cell  enclosing  the  point  (x,y)*  Determine  x*  ■  x-t, 
y*  “  r*j ,  the  (x,y)  location  with  respect  to  (l,j)«  Then  set: 

'\ 

T(x,y)  -  x'y*  fl+l>J+l  +  x'U-y’)  f1+l|j  ♦  (t-x*)  y*flfJ+1 

♦  (i-x»)  (1-y*)  ftj 
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Name:  CALCRES  (subroutine) 

Purpose:  Calculate  residuals  for  single  lcval  data* 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R*  Hoffman,  AER,  1986 

Referenced  by:  MASTOR2,  MAST0R6 

References:  None 

Commons  used:  DCONST,  DEL,  FGDATA,  UASIGXA  ■ 

Arguments:  PO 

POH  P£ 

KV  ^ 


I/O  units:  none 

Description: 

(0)  Residuals  are  returned  in  /DEL/* 

(1)  Set  t,  z  residuals  to  missing. 

(2)  Calculate  (u,v)  residuals.  If  (U^  ,  )  is  missing  set  (u,v)  residuals 

to  missing.  v  V 


CHLSKY 


Naas: 

CHLSKY  (subroutine) 

Purpose: 

Perforns  Cholesky  decoaposition  and  back  substitution  to 
solve  the  aatrix  equation  Ax“b  (Dahlquist  and  Bjorck,  1974, 
pp.  157-159;  Stobie,  1984). 

Author: 

D.  Norquist,  SASC,  1980  -  1986 

Docunentation : 

R.  Hoff nan,  AER,  1986 

Referenced  by: 

ASAP2 

References : 

none 

Coaaon8  used: 

none 

Argunents : 

FA  aatrix  A 

FB  r.h.s.  b 

N  size  of  systea 

X  solution  vector 

IFAIL  -  0  for  noraal  return 

■  1  if  a  pivot  saaller  than  IE-7  is  found. 

I/O  units: 

none 

Description: 

(1)  Decomposition  phsse:  Find  L  such  that  A  ■  LLT. 

(2)  First  back  substitution:  Find  y  where  Ly  ■  b. 

(3)  Check  for  pivot  too  snail*  Note  that  check  is  independent  of  data  b. 

(4)  Second  back  substitution:  Find  x  where  LTx  "  y. 


CNTOBS 


Name:  CNTOBS  (subroutine) 

Purpose  Counts  and  prints  out  the  nusber  of  valid  residuals 

classified  by  level  and  by  type* 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R*  Hoffman,  AER,  1986 

Referenced  by:  ASAP1 

References:  none 

Commons  used:  DCONST,  RESID 

Arguments :  NUM  N 

INUM 

ICHK.  -  0 

-  1 

I/O  units:  OUTPUT 

Description: 

(1)  Initialize:  write  heading,  zero  counters* 

(2)  Loop  over  all  observations,  switching  type  index  (IT)  as  necessary. 

(3)  At  each  level,  increment  height  counter  if  z*  is  present* 

(4)  At  each  level.  Increment  wind  counter  if  (uk,vfc)  is  present.  End  of  loop 

(2). 

(5)  Print  results.  Return* 


total  number  of  observations 

index  of  last  observation  of  types  1,  4,  2,  6 

before  consistency  checks 

after  consistency  checks 
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Mmi 

Purpose: 

Author: 

Document at ion : 


CQCV  (subroutine) 

Trans lutes  quality  indicators  froa  GWE  level  II  data  sets  to 
"generic  quality  uarks".  (See  Section  5*2.) 

D.  Norqulst,  SASC,  1980  -  1986 

R.  Hoffman,  AER,  1986 


Referenced  by: 

MAST0R1, 

MASTOR2 

References: 

none 

Commons  used: 

none 

Arguments : 

IQZ 

s,  quality  nark/indicator 

IQT 

T,  quality  mark/indicator 

IQQ 

q,  quality  ssark/indlcator 

I<** 

(u,v),  quality  aark/lndicator 

IDSI 

-  (not  used) 

ALAI 

-  (not  used) 

ALON 

-  (not  used) 

P 

-  (not  used) 

1/0  units :  none 


Description: 

(0)  On  input  IQZ  contains  the  ONE  quality  indicator,  on  output  it  contains 
the  ASAP  quality  mark.  The  saae  holds  for  IQT,  IQQ  and  IQW.  Variable 
IQC  contains  the  translation  table  froa  quality  indicator  to  quality 
nark*  It  is  based  on  WHO  (1978),  Appendix  A,  Table  IV,  Codes  IH,  IV  and 
IL.  IQC  is  therefore  valid  for  dsi.  -  11,  12,  21-24  only. 

(1)  Perfora  translation  for  z.  First  separate  two  digits  of  input  code  into 
a  horisontal  and  a  vertical  quality  control  indicator.  Translate  each 
separately  and  set  the  quality  mark  equal  to  their  maximum.  Note  that  a 
value  of  99  indicates  missing  both  for  a  quality  mark  and  a  quality 
indicator. 
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(2)  Perfora  translation  for  T. 

(3)  Perfora  translation  for  q. 

(4)  Perfora  translation  for  (u,v). 

(5)  Return* 
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Miam  * 

ESOBER  (subroutine) 

Purpose : 

Returns  estimated  observation  error  for  a  given  variable  for 

a  given  del  at  a  given  pressure  in  a  given  month.  Values  for 
satellite  heights  are  from  Norquist  (1986b,  pp.  36-37);  all 
other  values  are  from  Dey  and  Morone  (1985,  Table  5). 

Author: 

D.  Norquist,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman,  AER,  1986 

Referenced  by: 

ASaP2 

References : 

none 

Commons  used: 

none 

Arguments : 

IDSI  dsi 

P  P,  pressure  of  the  observation 

NV  variable  index;  1*Z,  2-u,  3-v 

PM  7  ,  pressure  at  the  mandatory  levels 

HANDLVL  L,  the  number  of  mandatory  levels 

EOV  EQ,  estimated  observation  error 

IMON  month  of  the  analysis 

Description: 

Cl)  If  (4  <  IMON  <  9)  move  June  values  to  fourth  column  of  EOE  table, 
otherwise  leave  February  values*  For  type  2  data  (A1REPS):  for 
del  ■  21,22,  set  EQ  -  4*9  a/s;  for  del  “  23,24,  set  EQ  ■  6.0  m/s. 

(2)  For  type  6a  data,  dsi,  -  61  (SATWINDS):  if  p  >  400  mb,  set  EQ  -  4.2  m/s; 
If  p  <  400  mb,  set  EQ  -  7.5  m/s. 

(3)  For  type  1  data,  dsi  -  13,14,15  (TWOS,  dropwindsondes),  set  EQ  -  2.0  m/s. 
For  dsi  ■  16  (constant  level  balloon),  set  E0  ■  5.0  m/s. 
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(4)  For  all  others,  interpolate  linearly  in  in  p  in  EOE  table  which  has  four 
categories  (coluans)  for  each  aandatory  level  (rows).  If  p  <  SO  ab,  the 
50  ab  value  is  used.  If  p  >  1000  ab,  the  1000  ab  value  is  used.  The 
four  categories  are 


Coluan 

Usage 

dai 

NV 

1 

RA0B  heights 

11 

1 

2 

RAO 3  winds 

11 

2.3 

3 

PIBAL  winds 

12 

2,3 

4 

SAT  heights 

41 

1 
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Nut:  FG  (subroutine) 

Purpose t  Bilinearly  Interpolates  first  quest  fields  to  the  observation 

location*  Returns  results  in  /FGDATA/and  in  arguments  ZSTAR, 
PO,  POH. 

Author:  D*  Norquist,  SASC,  1980  -  1986 

Docuaentatlon:  R.  Hoffman*  AER,  1986 


Referenced  by: 

MAS TORI 

,  MASTOR2,  MASTR04,  MAST0R6 

References: 

BIL1NR 

• 

Commons  used: 

CCONST , 

DCONST,  FGDATA,  FGFLDS,  SIGK 

Arguments : 

RLAT 

♦ 

RLON 

X 

ZSTAR 

z* 

PO 

aw 

Pk 

aw 

POH 

Pk 

IDS1 

dsl 

IWAT 

flag  for  R  computation  (aee  (3>). 

NRDLAT 

number  of  latitudes  which  have  been 

( X" "  !x  of  latitude  currently  stored  as  second  latitude  in 
/FGFLDS/*  That  is,  *2  -  -90  +  (NRDLAT  -  1)  A+«) 

I/O  units:  1,  5 

Description: 

(0)  Note  that  r  -  P* 

(1)  Determine  fine  grid  coordinates  (x,y)* 
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(2)  If  y  >  NRDLAT,  values  for  a  new  latitude  are  needed.  Shift  northern 
latitude  values  to  southern  latitude  values  in  /FGFLDS/.  Read  in  a 
new  northern  latitude,  incrementing  NRDLAT.  Go  to  (2). 

(3)  Interpolate  T,  U,  V,  R  from  /FGFLDS/  to  /FGDATA/ •  Relative  humidity 
la  Interpolated  only  if  data  is  type  1  (RAOBS)  or  if  data  is  type  4 
(SATEMS)  and  IWAT  *  0. 

(4)  Interpolate  P*  from  /FGFLDS/.  Set  Pk  and  P£  from  P*  and  /SIGK/. 

(5)  If  not  type  6  (i.e.  for  RAOBS,  AIREPS  or  SATEMS)  Interpolate  Z*  from 
/FGFLDS/  and  obtain  heights  by  integrating  T^  hydrostatically. 
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Nine: 


FLAGS  (subroutine) 


Purpose: 


Author: 


Perfora  gross  end  buddy  cheeks  on  the  data  stored  in  /RESID/ 
(see  Norquist,  1986b,  pp.  5-8).  Ail  comparisons  ere  made 
within  boxes  epproxiastely  10*  square  (see  Section  5.3). 
Gross  error  check  limit  is  ±3  forecast  error  standard 
deviations*  Buddy  check  follows  Bergman's  (1978,  1979) 
procedure  in  spirit,  with  the  addition  of  keep  flags.  Data 
which  fall  the  checks  are  set  to  missing. 

D.  Norquist,  SASC,  1980  -  1986 


Documentation:  R.  Hoffman,  AER,  1986 


Referenced  by: 

ASAP1 

References: 

none 

Commons  used: 

AS1AS2, 

Arguments : 

IBOX 

JDS  I 

QLAT 

QLON 

PTS 

NWO 

NUM 

CCONST,  DCONST,  IFLG, 

list  of  Bn 
list  of  d8in 
list  of 
list  of  Xn 

-  not  used 

-  not  used 
N 


RESID, 


SIGK 


I/O  units:  output 


Description: 

(0)  The  arrays  used  to  store  the  data  flags  and  indices  are  confusing.  In 

the  buddy  check,  within  a  box  at  a  single  layer  or  level,  observations  of 
the  same  type  are  compared  by  pairs  1,  j,  with  i  <  j.  I  will  use  the 
notation  in  this  documentation  that  F*  means  "i  puts  a  flag  on  j",  since 
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1  is  on.  top  of  j  in  the  symbol.  Also  I  use  WF  for  wind  flags  end  *F  for 
height  flags*  If  the  data  i  and  j  are  consistent,  then  ■  F^  »  0*  If 
they  are  inconsistent  an '  1  is  cf  better  quality,  then  F*  ■  1  and  Fj  •  0; 
if  they  are  of  equal  quality  then  F*  »  f|  »  1  and  If  j  Is  of  better 
quality  then  F^  -  1  and  p|  •  0* 

In  the  code  the  following  syabols  are  used 

N  -  1  IN  -  1  IFWI(IN,IN0)  -  IFZTI  (IN.INO)  -  *fJ 

NO  -  j  INO  -  j-i  IFWJ(IN,IN0)  -  WF*  IFZTJ  (IN.INO)  -  *pj 

(1)  Prediction  error  standard  deviations  Ep(Z4),  Ep(T4),  ^(U^)  and  Ep(Vt) 
for  the  mandatory  levels  in  5  latitude  belts  are  defined  from  Table  1  of 
Norquist  (1986b).  The  constants  a  and  b  for  the  buddy  check  (cf*  Eq.  7*1 
of  Bergman  (1979))  are  set  to.  3.0  and  1.5,  one  half  of  Bergman's  values. 

(2)  Begin  loop  over  the  18  10°  latitude  bands. 

(3)  Begin  loop  over  the  boxes  in  this  band. 

(A)  Calculate  the  box  index  for  the  current  box  and  begin  loop  over  all 
observations  in  the  box. 

(5)  If  array  overflow  occurs  (current  dimensional  constraint  is  100 
observations)  output  message  and  go  to  (6).  Otherwise,  add  current 
observation  to  the  list  of  observations  in  this  box.  NKB(i)  -  n(i),  i«l, 
Nfi  is  the  pointer  to  the  location  in  /RESID/  of  the  ith  observation  in 
the  current  box.  End  of  loop  (A). 

(6)  Loop  over  all  observations  in  the  box  1-1,  Ng,  extracting  residuals  and 
quality  marks  from  /RESID/  and  storing  them  in  local  arrays.  If  |+B|  > 
70s,  winds  are  converted  to  polar  stereographic  winds. 

(7)  Loop  over  layers  selecting  the  appropriate  kg  for  this  layer.  (See 
Dey  and  Morons,  1985,  p.  308).  If  ok  >  0.125  kh  -  2E-12a“2,  if 
0.125  >  ok  >  0.070,  then  kg  -  l.5E-12a”2,  and  if  0.070  >  ok,  then 
kh  -  1.0  E-12m~2. 
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(8)  Initialise  all  toaa  flags  and  keep  count tra  to  aero. 

(9)  Bogin  loop  over  all  obaarvatlona  la  tha  currant  box*  Determine  indices 
of  tha  aandetory  levels  above  and  below  tha  layer  and  level  pressures  of 
the  current  data*  That  is  find  tb»  it,  l*,  t*  such  that 

*  **  °^+1  >  **« 

Also  determine  the  latitude  belt  for  looking  up  the  forecast  error 
standard  deviations* 

(10)  Interpolate  Ep(Z£)  to  o£+j  and  Bp(Ug),  Ep(Vt)  to  linear  In  tn  p. 

Below  1000  nb  extrapolation  Is  used,  above  SO  ab  the  value  at  50  mb  is 
used. 

(11)  Perform  gross  error  check  on  current  observation  at  current  layer/level, 
keelduals  greater  than  3  forecast  standard  deviations  are  rejected ,  by 
changing  their  values  and  associated  quality  narks  to  missing*  A  list  of 
observations  In  the  current  box  which  have  failed  one  or  more  checks 
(gross  or  buddy)  is  maintained;  NCH(m)  for  a*l,  ICH  Is  the  ath_ 
observation  which  has  failed  a  check*  (See  (33)  and  (38).)  Note  that  at 
statement  label  72  NCHT  Is  1  If  a  height  error  has  been  detected  and  2  If 
a  wind  error  has  been  detected*  End  of  loop  (9)  on  observations. 

(12)  If  there  Is  more  than  one  observation  start  a  loop  over  all  observations 
1“1,  Ng.  This  loop  sets  the  Initial  flags*  Set  n(l)  *  NKB(i);  this  Is 
the  location  of  the  1th  observation  in  /RES1D/. 

(13)  Determine  the  coordinates  of  the  1th  observation  for  distance 
calculation.  A  polar  stereograyhic  projection  is  used  if  |+B|  >  70. 

(14)  Begin  loop  over  j«  1+1,  Ng.  Set  n(j)  -  NKB(j);  this  is  the  location  of 
Jth  observation  In  /RES1D/. 
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(IS)  Calculate  coordinate#  of  jth  observation  as  in  (13)  and  d^j ,  the  distance 
between  the  two  observations  as  in  ASAP1  (21-22)* 

Calculate  the  horlsontal  ss  correlation  (Day  and  Morone,  1985,  Bq*  2*8) 
and  the  vertical  correlation  (ibid.,  Eq.  2*9)  using  layer  pressures* 

If  or  Zj  is  Biasing  go  to  (19)*  Otherwise  calculate  the  vertical 
correlation  (Eq*  2*9)  using  level  pressures,  and  the  value  of  p**  of  the 
zz  total  correlation  (ibid.,  Eq.  27)..  Calculate  the  r.h.s.  and  i.h.s.  of 
the  buddy  check  criterion  for  a  (Bergman,  1979,  Eq.  7.1):  these  are 
(a-bp**)Ep(x)  and 

Determine  QL(z1)  and  (^(Zj),  the  z  quality  levels.  (See  Section  5*2.) 

If  the  buddy  check  criterion  is  violated  for  c,  then  set  toss  flags: 

(a)  If  QL(zi)  <  QL(*j)  then  -  1  and 

(b)  If  QL(Zj)  <  QjUj)  then  *F^  -  1. 

If  the  criterion  is  not  violated  and  if  | p**  |  >  0.75,  then  increment  the 

keep  counters  K^,  Kj. 

(19)  If  either  wind  is  missing  go  to  (22).  Calculate  |p“j|  and  |p]^|  using 
equations  from  Dey  and  Morone  (1985;  either  (A6)  and  (A9)  or  (B6)  and 
(B9)  as  appropriate). 

(20)  Determine  wind  quality  levels  QJj(u1,v1),  Q^u^Vj).  Q^Cuj.Vj), 

QL(uj,vj). 

(21)  Set  flags  for  winds  as  in  (18).  The  buddy  check  criterion  is  violated  if 
either  u  or  v  components  violate  their  separate  tests.  The  keep  counters 
are  incremented  only  if  |p^j|  and  |p^j  ars  both  >  0.75. 

(22)  End  loop  (14)  on  j.  End  loop  (12)  on  1. 


(16) 


(17) 

(18) 
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(23)  Loop  over  all  observations  in  tha  box*  Remove  all  toaa  flags  on 
observations  which  hava  2  or  more  kaap  flaga*  l.a.  If  (*1^  >  1)  chan 
*f|  -  0  for  J  "1*  N|  and  similarly  for  tha  winds* 

(24)  Initialise  Interatlva  rajacdon  cycle*  The  alniaui  number  of  toaa  flaga 
which  require  an  observation  to  be  toaaed  (Nj)  la  aet  to  2*  unlaaa  there 
are  only  2  observations  In  a  box  in  which  caae  It  la  aet  to  1* 

(25)  Begin  rejection  cycle*  Initialise  sums  (W$A,  8SiP  NR).  Begin  loop  over 
all  obaervationa  In  box*  i-l*Ng*  (Arraya  1QZ  and  IQW  are  used  here  for 
wSt  and  ‘Sj*) 

Accuaulate  the  number  of  toaa  flaga  on  1  for  j  >  1*  (Process  f|  for 
j  >  1). 

Accuaulate  the  number  of  toss  flaga  on  j  by  1  for  j  >  1*  (Process  Fj  for 
j  >  1).  Note  (25-27)  are  equivalent  to  creating  the  sums 

■»  . 

S*  ■  \  F~  for  all  1  and  for  s  and  (u*v). 
j-1 

If  WS^  <Nj  and  *S1  <  Nj.  then  Increment  N^  the  number  of  observations 
which  are  assured  of  being  kept*  End  loop  (25)  on  1. 

(29)  II  **g  ■  Ng  then  all  observations  remaining  are  acceptable:  go  to  (35). 
Otherwise  there  are  more  to  be  tossed* 

(30)  Determine  minimum  and  maximum  number  of  toaa  flags  (aln  Sit  max  St)  and 
the  index  iMI  which  corresponds  to  the  observation  with  the  maximum,  for 
both  x  and  (u*v),  separately* 

(31)  If  max  Sj:  in  >  4  truncate  to  a  multiple  of  4  (e*g.  7  becomes  4)  for  both 
x  and  (u,v),  separately. 


(26) 


(27) 


(28) 
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(32)  Toss  data  (by  setting  the  local  residuals  to  Biasing  )  for  those  1  which 
have  >  the  (truncated)  aax  S^.  Maintain  a  list  of  pointers  to  the 
observations  in  the  local  arrays  which  have  been  tossed  during  this 
rejection  cycle.,  e.g.  NZ(a),  a»l,  INTZ  contains  the  pointers  to  the 
tossed  z  observations. 

(33)  Update  the  list  of  pointers  to  observations  which  have  failed  one  or  aore 
tests  (NCH(n),  a-1,  1CH)  to  reflect  changes  aade  in  (32)  using  the  NZ  and 
N W  lists. 

(34)  Remove  all  toss  flags  iaposed  by  observations  which  have  just  been  tossed 
on  other  observations.  That  is  for  i-1,  Mg  if  >  (truncated)  aax  (S^) 
then  set  Fj  -  0  for  j«l,  NR,  for  z  and  (u,w)  separately. 

(35)  Go  to  (25)  to  start  next  iteration  unless  10  iterations  have  already  been 
performed. 

(36)  Loop  on  observations  in  the  box  i«l,Ng.  Restore  1QZ,  IQW  arrays  froa 
/RESID/.  These  arrays  (Q(z),Q(u,v))  were  used  to  store  and  WS*. 

Since  some  observations  have  been  tossed,  set  Q(z)  to  aissing  If  the  s 
residual  is  aissing,  and  set  Q(u,v)  to  aissing  if  the  (u,v)  residual  is 
missing.  Check  if  any  observations  still  aight  be  tossed.  If  this 
occurs,  print  message. 

(37)  End  of  loop  (7)  on  layers. 

(38)  Move  changes  back  to  OBD  array,  using  list  NCH.  If  |tB|  >  70  convert 
winds  from  polar  stereographic  back  to  (X,  $)  coordinates. 

(39)  End  of  loops  (2)  and  (3)  on  buddy  check  boxes.  Print  message  and  return. 
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Nmm: 


Purpose: 


Author: 


LOWTMP  (subroutine) 

Converts  analysed  residuals  *k  to  temperature  residuals  tk 
and  corrects  the  first  guess  value  Tk  .  (See  Norqulst, 
1986b,  p.  18,  pp.  55-58.) 

D.  Norqulst,  SASC,  1980  -  1986 


Documentation:  R.  Hoffman,  AER,  1986 


Referenced  by:  ASAP2 


References:  none 

i 

Commons  used:  DCONST,  LOUT 


Arguments : 


CORZ  *kf  layer  analysis 

CORZU  *u,  underground  layer  analysis 

TG  Tk  ,  first  guess  on  Input,  updated  analysis  on 

output 

SL  scaled  pressures  pk 

SLU  scaled  pressure  pu  in  the  layer  underground 

Note:  all  arguments  are  variables  for  the  grid  point  being 

analysed. 


I/O  units:  none 

v 

Description: 

(0)  Vertical  indexing  used  In  this  subroutine  is  illustrated  in  the  following 
figure.  The  first  value  is  used  in  (1)  and  (4).  The  second  value  is 
used  otherwise. 


LOWTHP 


k  "  K+1,K 


k  -  K,K-1 


top  of  the  atmosphere 


k  -  3,2 


k  -  2,1 


k  -  1 


surface 


(1)  Compute  T  corrections  between  two  layers  (i.e.,  at  "interfaces" ), 

-  Vi>/Un  pk  "  to  "k-i*  k  ' 2 . V1 

These  temperatures  must  be  multiplied  by  -g/R  to  become  diaensionalized 
(see  (3)). 

(2)  Reset  index  so  that  t*  -  t*  ,  k  ■  1,...,N, 

k  k+1  k 

(3)  Apply  Flattery  algorithm,  treating  tA  as  layer  temperatures*  Since  the 
sigma  structure  /S1GK/  doesn't  change,  the  required  constants  are  stored 
in  /LOWT/.  First  form  BT\  then  A*BTA  and  finally  (ATA)“1AT3Ta* 

Here  TA  is  the  vector  of  tA>  Dimensionalize  the  result  by  ~g/R  to  get  T, 
the  luyer  temperatures. 

(4)  Set  uppermost  layer  residual  to  zero:  tR+1  «  o.  (Note:  layers  include 

subsurface  layer.)  Update  temperatures  +  t^+i  k  *»  1,K. 
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Naas;  LOWTMPS  (subroutine) 

Purpose:  Implements  the  Flattery  algorithm  for  Interpolating  layer 

teaperaturea  to  level  temperatures  by  the  method  of  least 
squares.  This  version  assuaes  the  topmost  pressure  Is 
greater  than  ssro.  There  should  be  no  missing  values. 
(Sea  Norqulst,  1986b,  Appendix  A,  pp.  55-58.) 

Author:  D.  Norqulst,  SASC,  1980  -  1986 

Documentation :  R.  Hoffman,  AER,  1986 


Referenced  by: 

SATLTMP 

, 

References : 

MATRIX 

Commons  used: 

iDCONST 

Arguments : 

TLAY 

Tt 

ML 

L  -  1 

TL 

T*  ,  level  temperature 

PL 

P*  (in  practice  these  are  scaled.) 

ATA 

<!**>-» 

MLPI 

L 

1/0  Units: 

none 

Description: 

(0)  Note  that  the  aaxlmua  number  of  levels  Is  20. 

(1)  Set  scale  for  temperatures  to  -g/R.  Calculate  weights  for  interpolating 
layer  to  level  temperatures  for  levels  2,...,L-1.  Note  that  the  sum 
of  the  two  weights  is  one.  Layer  pressures  are  assumed  to  be  given 
by  P,  -  (P ;  PJtl)1/2  • 
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(2)  Form  A^  matrix. 

(3)  Form  (A^A)“*  matrix. 

(4)  Form  BTt  first  scaling  the  temperatures.  Note  special  handling  of  BT(L) 
and  BT(2L-1)  corresponding  to  levels  1  and  L. 

(5)  Form  AT(BT). 

(6)  Calculate  level  temperatures  (ATA)“*  (A^BT)  and  diaensionalice. 
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Nun:  MASTORl  (subroutine) 

Purpose:  Reeds  type  1  dsts  (RAOBS)  froe  unpecked  GWE  Level  II  data 

set,  interpolates  to  sigee  end  calculates  residuals,  which 
are  then  stored  in  /RES ID/.  (See  Norquist,  1986b,  pp.  31- 
33.) 


Author: 

Docuaentation: 

Referenced  by: 
References: 


D.  Norquist,  SASC,  1980  -  1986 
R.  Hoff aant  AER,  1986 
AS  API 

CQCV,  FC,  PTOSIG 


Coaeons  used:  DCONST,  FGDATA,  RESID,  U AD AT A,  UASIGMA 


Arguaents :  IDATE 

ITIME 
NV 
HUM 
QLAT 
QLON 
JDSI 
1BOX 

I/O  Units:  2, OUTPUT 


date 

time 

-  (not  used ) 

N 

list  of  *n 
list  of  xn 
list  of  dsln 
list  of  Bn 


Description: 


(1)  Initialise  tiee  and  date  window  as  A  3  hours  about  date/tine  which  was 
input.  Set  ICHNG  ■  0;  this  is  a  flag  used  to  overwrite  a  previous 
observation  with  a  new  observation  of  higher  quality.  Initialize  N,  L. 


(2)  Read  a  record  froa  unit  2,  using  format  described  by  Norquist,  (1984, 
pp.  18-19).  On  EOF  go  to  (19). 
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(3)  If  not  acceptable  go  to  (2).  Specifically  dsl  ■  75,  cloud  data,  are 
skipped.  The  presence  of  other  non-type  1  data  results  in  an  ABEND. 
Records  with  missing  p,  negative  z  or  date/time  not  in  Che  window  are 
skipped. 

(4)  If  X,  4  or  dsi  do  not  match  xN,  +N,  or  dsiH>  go  to  (12).  That  is  if 

current  record  is  not  part  of  current  observation,  it  is  time  to 

complete  processing  the  current  observation  and  begin  new  one.  Note 
that  the  test  makes  use  of  fact  that  X,  +  are  coded  only  to  hundreths  of 
degree. 

(5)  If  level  dimension  (currently  90)  is  exceeded  in  /UADATA/ ,  skip  this 
level,  go  to  (2). 

(6)  If  this  is  not  the  first  level  for  the  current  observation,  i.e..  If 

l  *  1,  go  to  (8).  Set  XN,  4n,  dslfl  to  X,  4,  dsi.  Calculate  BQ  (see 

Section  5.3). 

(7)  If  N  >  1,  then  loop  over  all  previous  observations  searching  for  a 

matching  Xn,  If  match  is  found  then:  If  dsi  <  dsl^,  then  the  new 

observation  is  better  (see  UNO  (1978),  Appendix  A,  Table  1).  Set 
ICHNG  -  1  to  force  overwriting  observation  n.  (To  do  this  temporarily 
set  N  -  n.)  Go  to  (8).  Else  if  dsi  >  dain»  then  flush  current  records 
until  new  location  is  found.  Go  to  (3). 

(8)  Translate  T  in  Celsius  to  T  in  Kelvin,  T  to  Q,  wind  speed  and 

d 

direction  to  (U,V).  If  quality  indicator  is  99  for  any  variable,  set 
corresponding  value  to  missing  (DNN  -  -999.9).  Winds  at  the  pole  are 
set  to  missing. 

(9)  Combine  quality  indicators  and  translate  into  qualtiy  marks  (CQCV). 

(10)  Set  all  quality  marks  >  1  to  99,  set  corresponding  values  to  missing. 
This  keeps  only  data  which  have  been  checked  satisfactorily  and  data 
which  have  not  been  checked. 
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(11)  Increment  L.  Go  co  (2). 

(12)  Save  X*  4  for  Che  new  observation  which  has  been  detected*  If  L  ■  1 
(which  can  happen  only  if  N  -  1),  go  to  (6)* 

(12a)  Eliminate  duplicate  and/or  out  of  order  levels  from  current  observation. 

(13)  Calculate  first  guess  values  at  XN*  for  ZA,  1*^, 

and  (FG).  Although  IWAT  -  0,  is  calculated*  since  this  is  type  I 
data* 

(14)  If  the  observation  has  only  a  single  level- skip  it  by  going  to  (17a). 
Interpolate  observation  to  sigma  coordinates  (PTOSIG).  If  all  sigma 
coordinate  data  its  missing  go  to  (17a).  Convert  to 
If  Qk  or  is  missing*  then  is  missing. 

(15)  Compute  residuals.  Make  sure  that  missing  data  results  in  missing 
residuals. 

(16)  Store  residuals  and  associated  quality  marks  in  /RESID/.  (At  this  point 
quality  marks  are  either  0  for  checked  residuals*  1  for  unchecked 
residuals*  or  99  for  missing  residuals.) 

(17)  Restore  N*  ICHNG  if  overwriting.  Increment  N  for  next  observation. 

(17a)  If  EOF  go  to  (19).  Otherwise  more  contents  of  current  record  to  start 
of  /UADATA/  arrays*  since  this  data  is  start  of  the  new  current 
observation.  This  resets  L  •  1.  Go  to  (6). 

(18)  Position  unit  2  at  start  of  next  file.  Go  to  (12a) 

(19)  Undo  (17):  N  ■  N-l.  Return. 
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Name : 

MAST0R2 

(subroutine) 

i 

Purpose: 

Reads  type  2  data  (AIREPS,  etc*)  from  unpacked  WE  Level  II 

data  set 

.  From  each  report  (u,v)  residuals  are  calculated 

i 

] 

( 

i 

and  assigned  to  the  analysis  layer  doseat  to  the  reported 

height* 

(See  Norquist,  1986b,  pp*  33-34*) 

] 

j 

1 

Author: 

D.  Norquist,  SASC,  1980  -  1986 

j 

Documentation: 

R.  Hoffman*  AER,  1986 

] 

Referenced  by: 

ASAP1 

1 

1 

1 

References: 

CALCRES, 

CQCV,  PC 

j 

1 

Commons  used: 

DCONST, 

DEL,  FGDATA,  RESID,  UADATA,  UASIGMA 

.  1 

J 

j 

j 

Arguments : 

IDATE 

date 

i 

i 

ITIME 

tine 

NW 

-  (not  used) 

NUM 

N 

) 

i 

\ 

QLAT 

list  of  #n 

QLON 

list  of  Xn 

JDS  I 

list  of  dsin 

\ 

I  BOX 

list  of  Bn 

\ 

| 

I/O  Units: 

2 

I 

Description: 

/ 

(1)  Initialize  date/tine  window  aa  *  3  hours  about  input  date/time. 

Initialize  N  for  next  observation*  Initialize  /UASIGMA/  humidities 
to  missing  (Q^)  .  Initialize  all  of  /UADATA/  to  missing  (P^,  ,  U^, 

Vf  V  Z£)#  Initiali*e  a11  of  /FGDATA/  to  missing  (Z£,  TR,  Uk,  Vk,  Rj,)  . 


! 
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(2)  Read  a  type  2  header  record.  See  Norquist,  (198A,  pp.  20-21).  Ac  EOF 
go  Co  (23). 

(3)  Check  if  data  is  acceptable*  Non-type  2  data  results  in  an  ABEND.  Data 
outside  date/tlme  window  are  skipped  by  going  to  (22). 

(4)  Set  Xjj,  and  dsl^  to  values  read*  Calculate  (see  Section  5.3). 

(5)  If  current  observation  is  collocated  with  any  previous  observation, 
including  those  of  other  data  types,  the  current  observation  is  skipped 
by  going  to  (22). 

(6)  Store  Z  and  T  and  associated  quality  narks  from  the  header  record, 
converting  degrees  Celsius  to  degrees  Kelvin. 

(7)  Loop  on  wind  records  associated  with  current  header  record  which  was 
read  in  (2).  (Currently  sorted  data  always  have  a  single  wind  per 
header  record*)  First  wind  report  is  preset  to  missing,  then  data  are 
skipped  if  quality  nark  indicates  wind  information  is  missing  or  if 
location  is  at  the  pole*  Otherwise  wind  speed  and  direction  are 
converted  to  (U,  V)  and  the  quality  mark  is  saved* 

(8)  Skip  the  optional  record  if  present* 

(9)  Translate  the  quality  indicators  to  quality  marks  (CQCV). 

(10)  Only  quality  marks  of  0  and  1  are  retained*  If  Q(Z)  >  1,  skip  this 
observation  by  going  to  (2)*  If  Q(T)  >  1  or  Q(U,V)  >  1,  the  associated 
values  and  quality  narks  are  set  to  missing. 

(11)  Interpolate  first  quess  to  XN,  +N  (FG).  Note  that  R  is  not 
interpolated* 

a# 

(12)  If  Z  <  ZA  or  Z  >  Zg,  then  go  to  (21).  Otherwise  determine  the  Index  kc 
of  the  Z£  closest  to  Z  • 
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(13)  Calculate  the  pressures  associated  with  zr  and  Z  using  formulas  for  a 
standard  atmosphere  as  described  by  Norqulst  (1986b,  Eq.  24  and 

0^0  as 

subsequent  discussion).  Call  these  P(2)  and  P(Z). 


(14)  The  actual  pressure  associated  with  the  observation  la  then  taken  to  be 
p  defined  so  that  p  -  P£  +j  *  P(2)  -  P(Z)  • 

(Note  that  i*  is  at  pressure  p£+^.) 

(15)  Determine  the  sigma  layer  kv  which  has  pressure  closest  to  p. 

(16)  If  T  is  present,  interpolate  first  guess  linear  in  in  p  to  p  and  then 

add  residual  to  Tk  to  obtain  Tk  •  This  effectively  assigns  the 
%  *v 

residual  to  layer  ky.  if  T  is  not  present  'set  T|^  to  missing. 

(17)  If  (U,V)  is  present  interpolate  first  guess  linear  in  in  p  to  p  and  then 
add  the  residual  to  (U^  ,Vk  )  to  obtain  (u^  ,Vj^).  This  effectively 
assigns  the  residual  to  layer  ky.  if  (u,V)  is  not  present, 

set  (u.  ,v.  )  to  missing. 

(18)  For  all  k  *  ky,  set  T^,  Uk,  Vk  to  missing. 

(19)  Compute  residuals  (CALCRSS). 


(20) 


Store  results  in  /RESID/ • 
Increment  N. 


Each  observation  contains  only  P  ,  u^  ,  • 


(21)  If  only  a  single  wind  is  associated  with  the  current  header  go  to  (2). 
This  will  always  be  the  case  with  sorted  data.  Otherwise  process  other 
wind  data  roughly  following  steps  (4),  (11),  (12),  (13),  (14),  (15), 
(12),  (19),  (20).  Then  go  to  (2). 

(22)  Skip  wind  record(s)  and  the  optional  record  if  it  is  present,  which  are 
associated  with  current  header  record. 


(23)  On  EOF  position  unit  2  after  the  EOF  mark.  Undo  last  N  increment  and 
return. 
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Name: 


MAST0R4  (subroutine) 


Purpose:  Reads  In  SAXTRMS  fros  GWK  level  II  unpacked  data  set  which 

doss  not  Include  SATTRMS  over  land*  Data  Is  transformed  to 
sigma  coordinate  height  residuals  and  stored  in  /RESID/. 
(See  Norqulst  1986b,  pp*  34-36*) 


Author:  D.  Norqulst,  SASC,  1980  -  1986 

Documentation:  R*  Hoffman,  AER,  1986 


Referenced  by:  AS API 


References:  FG,  SAILTMP 


Commons  used:  D00NST,  FGDATA,  RESID,  U AD AT  A,  UASIGHA 


Arguments: 


IDATE 

date 

ITIME 

time 

NW 

- 

NUM 

N 

QLAT 

list 

of 

^n 

QLON 

list 

of 

Xn 

JDS  I 

list 

of 

daln 

IBOX 

list 

of 

Bn 

(not  used) 


I/O  Units:  2 


Description: 


(1)  Initialise  date/time  window  as  ±  3  hours  about  date/time  as  input* 

Increment  N  for  next  observation*  Initialise  observations  interpolated 
to  sigma  as  missing  for  U,  V  and  Q,  i.e.,  Vk«  Qk-  DNN.  Initialise 
all  first  guess  values  at  observation  location  to  missing;  i.e*, 

«W  M  W  Ml 

K  -  V  V  V  Tk  - 
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(2)  Read  in  type  4  data  record.  This  should  be  the  group  header  record  for  a 
number  of  layer  records  (see  Norquist,  1984,  pp.  25  -  27).  On  EOF  go  to 
(18) 

(3)  If  not  acceptable  go  to  rejection  loop.  Specifically  dal  must  be  41  and 
the  date/tiae  must  be  acceptable.  Note:  soundings  over  land  are  not 
included  in  the  input  data  set. 

(4)  Set  XN,  +n,  dsiN  to  values  read.  Calculate  Bg* 

(5)  If  collocated  with  a  previous  observation,  possibly  of  another  type,  go 
to  rejection  loop  (17). 

(6)  Set  quality  mark  to  0  if  IQCI  (WHO  (1978),  Appendix  A,  Table  38)  and  ISCR 
(WHO  (1978),  Appendix  A,  Table  37)  indicate  obaervatlon  was  found  correct 
during  quality  control.  Quality  nark  is  set  to  1  if  observation  was  not 
checked.  All  other  observations  are  skipped  by  going  to  rejection  loop 
(17).  Note  that  this  specifically  excludes  observations  corrected  during 
Special  Effort  activities. 

(7)  Begin  reading  layer  data  records.  See  Norquist  (1984,  pp.  26-27)  for 
details  on  the  four  types  of  layer  inforaatlon.  MAST0R4  assumes  that 
thickness  and  mean  temperatures  are  not  both  present  and  that  standard 
and  nonstandard  layer  precipltable  water  are  not  both  present.  Obtain 
first  guess  values  for  Z*,  P^,  P£,  T^,  U^,  Vk,  Z£,  and  if  layer 
precipltable  water  is  present,  R^  • 

(8)  Xf  no  thickness  data  is  present  go  to  (12).  Read  data  into  arrays  of  P^ 
and  Z.  where  P#  is  the  pressure  at  the  top  of  layer  i,  z.  1®  the  height 
at  relative  to  Pj.  Thus  Zj  -  0  •  The  layer  temperatures  T^  are  then 
calculated  hydrostatically  (see  Section  5.1.2).  Go  to  (10). 

(9)  For  mean  temperature  data  ?t  «nd  «re  read  in.  Data  associated  with 
missing  pressures  or  with  layers  of  zero  or  negative  mass  (Ap)  are 
skipped.  The  pressure  at  the  lower  boundary  should  match  the  previous 
upper  boundary  pressure.  If  it  is  too  large  an  error  exists,  if  too 
small,  then  a  layer  with  missing  is  defined  between  the  two 
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prc««urM<  Indexing  used  her*  In  normal:  T  la  the  Man  temperature 
between  and  *A+j» 

(10)  Interpolate  to  (at  algna  layers)  and  Z£  (at  algae  levels)  using 
Flattery  algorithm  and  hydrostatic  relationship  along  with  some  first 
quaes  information  (SATLTMP).  Mote  dlMnslonal  limit  of  20  layers  in 
matrix  A 1a. 

(11)  Set  quality  Mrks  for  Z£  and  equal  to  the  quality  Mrk  defined  in  (6) 
for  all  levels  which  have  data  and  to  99  for  missing  data*  The  layer  by 
layer  quality  marks  of  the  Level  II  data  are  ignored*  Go  to  (13)* 

(12)  Thickness  data  Is  not  present;  set  all  quality  Mrks  and  all  T^,  and 
Z£  deta  values,  to  missing  in  /UASICMA/* 

(13)  Skip  records  corresponding  to  layer  predpitable  water  data*  Set  all 
associated  quality  Mrks  and  data  values  to  aiaslng  in  /UASIGMA/ •  If 
Man  temperature  data  is  to  be  read  next,  then  go  to  (9). 

(14)  Begin  computation  of  residuals*  Store  P*  in  /RESIT)/.  Code  is  present 
here  to  convert  Q  to  R  (see  Section  5.1.4)*  However,  Q  will  always  be 
missing  (see  (13))* 

(15)  Calculate  residuals;  if  data  is  missing  In  /UASIGMA/,  the  quality  marks 
and  residuals  are  set  to  missing.  Hind  quality  marks  and  residuals  are 
autoMtically  set  to  missing*  All  residuals  are  stored  in  /RESID/. 

(16)  IncreMnt  N  and  go  to  (2)* 

(17)  Rejection  loop:  Read  thru  all  records  of  the  current  sounding. 

Go  to  (2).  Unexpected  EOF  results  in  an  ABEND  here. 

(18)  EOF  handling  (see  (2)). 
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Name:  MAST0R6  (subroutine) 

Purpose:  Reads  in  type  6a  data  (SATWINDS)  f roa  unpacked  GUI  Level  11 

data  set.  Froa  each  report,  (u,v)  residuals  are  calculated 
and  assigned  to  the  signs  layer  closest  to  the  reported 
pressure.  (See  Norquist,  1986b,  pp.  41-42.) 


Author: 

0.  Norquist 

,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman, 

AER,  1986 

Referenced  by: 

ASAP1 

References : 

CAJ.CRES,  FG 

Commons  used: 

DCONST,  DEL 

,  FGDATA,  RESID,  UADATA,  UASIGMA 

Arguments : 

IDATE 

date 

ITIME 

time 

NH 

-  (not  used) 

NUM 

N 

QLAT 

list  of  dn 

QLON 

list  Of  Xy 

JDSI 

list  of  d8in 

IBOX 

list  of  Bn 

I/O  Units: 

2,  OUTPUT 

Description: 

(1)  Set  up  quality  nark  translation  tables  1CP  and  ICW  (see  (7)).  Initialize 
date/time  window  as  ±  3  hours  about  input  date/tiae.  Increaent  N  for 
next  observation. 

(2)  Initialize  /UASIGMA/  humidities  to  missing  (Q^)  .  Initialize  all 

/UADATA/  to  missing  (Pt,  Tt,  U^,  V^,  Qt,  Z t)  .  Initialize  all  /FGDAIA/ 
t  o  mi  s  s  i  ng  (  ,  T^ ,  U^ ,  )  . 
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<3)  Read  a  type  6a  record.  See  Norqulst  (1984,  pp.  29-30). 

(4)  Check  if  data  is  acceptable.  Non-type  6a  data  results  In  an  ABEND*  Data 
outside  date/tiae  window  are  skipped  by  going  to  (3). 

(5)  Set  dN,  \y,  daiM  to  values  read.  Calculate  (see  Section  5*3), 

(6)  If  current  observation  is  collocated  with  any  previous  observation, 
including  an  observation  of  another  data  type',  the  correct  observation  is 
skipped  by  going  to  (3).  Note  collocation  criterion  agreement  to  0*1 
degrees  in  latitude  and  longitude  since  for  these  data  location  is  coded 
to  0.1  degrees,,  but  other  data  are  coded  to  0.01  degrees* 

(7)  Sec  quality  marks  Q(T)  and  Q(U,V)  from  quality  indicators  IQCP  and  1QC 
read  at  (3)  using  the  translation  tables  ICP  and  ICW.  These  tables 
translate  VHP  (1978)  Appendix  A  Tables  30  and  36  into  the  ASAP  "generic 
quality  marks"  (see  Section  5*2).  Note  that  values  of  IQCP  or  IQC  >  3 
are  translated  to  a  quality  mark  of  5.  However  Table  36  has  been 
extended  by  the  USA  Special  Effort  (U*  of  Wisconsin)  to  Include 
indicators  ranging  from  5-9.  Data  with  IQC  indicators  6  and  7  should  be 
translated  to  a  quality  nark  of  0,  5  to  1,  8  to  6  and  9  to  8 „ 

(8)  If  T  is  minding  or  Q(T)  >  1,  set  T  and  Q(T)  to  missing.  Otherwise 
convert  degrees  Celsius  to  degrees  Kelvin. 

(9)  If  wind  speed  or  direction  is  missing  or  Q(u,v)  >  1,  set  (u,v)  and 
Q(u,v)  to  missing.  Otherwise  convert  wind  speed  end  direction  to  (u,v). 

(10)  If  T  and  («,v)  are  missing  skip  this  observation  by  going  to  (3)* 

(11)  Interpolate  first  quess  to  XN,  $N  (FG).  Note  that  R  is  not  interpolated. 

(12)  If  reported  pressure  P  is  not  between  top  and  bottom  sigma  layer 

•v  <v 

pressures,  (Pj,P^)  then  skip  the  current  observation  by  going  to  (3). 
Otherwise  determine  ky,  the  index  of  the  sigma  layer  pressure  closest  to 
the  observation  pressure* 
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(13)  tf  T  is  present,  interpolate  first  guess  linear  in  in  p  to  F  and  then  add 
residual  to  Tk  to  obtain  T.  .  This  effectively  assigns  the  residual  to 
ky.  If  T  is  oissing,  set  Tk  to  missing* 


(14)  If  (u,v)  is  present  interpolate  first  guess  linear  in  in  p  to  F  and  then 
add  the  residual  to  (U.  ,  V.  )  to  obtain  (U^  ,  V.  )  .  This  effectively 
assigns  the  residual  to  ky.  If  (u,v)  is  Biasing,  set  (l^  ,Vk  )  to 
missing. 


(15)  For  all  k  ♦  ky,  set  T^,  u^,  v^  to  missing. 


(16)  Compute  residuals  (CALCRES). 


(17)  Store  results  in  /RESID/.  Each  observation  contains  only  F  ,  u^  ,  vk 


(18)  Increment  N  and  go  to  (3)  unless  N  now  exceeds  6000,  the  current  array 
dimension  in  /RESID/. 


(19)  EOF  or  dimensions  exceeded.  Close  unit  2.  Undo  last  N  increment. 
Return. 
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N«ae: 

POINTS  (subroutine) 

Purpose : 

Calculates  the  number  of  data  quantity  points  (?TS)  for  each 
observation*  Data  quulity  of  SAtEMS  is  Accounted  for. 
Additional  credit  is  given  If  z  and  (u,v)  are  both  prerent. 
(See  Horqulst,  1986b,  pp.  8-9.) 

Author : 

D.  Norquist,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman,  AER,  1986 

Referenced  by: 

AS  API 

References: 

none 

Commons  used: 

DCONST,  RESID 

Arguments : 

JDS1  dSln 

PTS  polntsQ 

NUM  N 

I/O  units: 

none 

Description: 

(1) 

Loop  over  all  observations  n»l,  N.  Initialize  polntsn  -  0. 

(2) 

Set  factor  to  1.0  unless  dsi-/10*4  (i.e.  SATEMS)  in  which 

case  factor  is  set  to  the  ratio  of  the  observational  standard 

deviations  of  RAOBS  to  SATEMS  (0.42  in  winter  and  0.53  in 

suemer). 

(3) 

Loop  over  all  layers  k  ■  1,  K. 

If  z£n  Is  present  add  factor  to  pointsn. 

If  (u,v)fc„  is  present  add  factor  to  points.,* 

If  *£n  and  (u,v)i,^  are  present  add  factor  to  points,,. 

(A) 

End  loop  (3)  on  k.  End  loop  (1)  on  n.  Return. 

67 


PTOSIG 

Name:  PTOSIG  (subroutine) 

Purpose:  Interpolates  RAOBS  to  signs*  Assumes  that  tenperature  varies 

linearly  in  in  p  and  hence  that  height  varies  quadratlcally 
in  In  p.  Both  T  and  Z  observations  sra  used  to 

m*  Mi 

calculate  Tfc,  the  signs  layer  tenperatqres*  The  t^ 
supplemented  by  first  guess  heights  as  needed  are  used  to 
calculate  the  Z*  .  Hinds  and  in  (specific  tumidity)  are 
interpolated  linearly  in  inp.  (See  Norquist  1986b,  pp.  32- 
33,  61-63.) 

Author:  0*  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER.  1986 

Referenced  by:  MAST0R1 

References:  none 

Commons  used:  DCONST,  FGDATA,  UADATA,  UASIGHA 

Arguments :  NLVL  L 

ZSTAR  Z* 

\ 

POH  Pk 

I/O  units;  none 

Description: 

(1)  Begin  loop  over  all  sigma  levels  k-1,  K+l .  This  loop  interpolates  Z£  and 

T^  to  a  preliminary  estimate  of  Zr  •  Find  £p  the  first  1“1,  L  such 
that  P^  <  .  Go  to  (3). 

(2)  Fall  borough  loop*  If  P^  is  within  0.1  nb  of  PK+j,  allow  extrapolation 
of  z  by  setting  -  L.  Otherwise  go  to  (8). 
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(3)  If  •  1  go  Co  (8).  (No  downward  extrapolation.) 

(4)  If  (Pti-l  -  Pt)  >  300,  go  to  (8). 

(5)  Find  ft^,  the  first  level  ft  at  or  below  ft^-l  which  has 

both  Z^  and  data.  Find  ftu,  the  first  level  ft  at  or  above  ft^  which  has 

both  Z&  and  data.  If  either  can't  be  found,  go  to  (8). 

(6)  If  (P,  -  P.  )  >  300  mb,  go  to  (8). 

Xft  *u 

(7)  Interpolate  Z  and  T  quadratically  in  ftn  p  to  obtain  the  preliminary 

*  *  mm 

estimate  of  Z^  .  Note  that  here  is  the  height  at  P^,  not  at  Pk+j  • 
This  algorithm  is  comprised  of  Eqs.  C7,  C3,  C6,  C2  and  Cl  of  Norqulst 
(1986b,  pp.  61-62),  except  that  the  constants  A  and  B  are  multipled  by 
(-1)  everywhere.  The  quality  mark  Q(Z£)  is  set  equal  tu  the  maximum 
of  the  quality  marks  associated  with  the  four  input  data, 

T£  and  Z^  at  ft&  and  ftu>  Go  to  (9). 

(8)  Do  not  interpolate.  Set  Z^  to  missing. 

(9)  End  of  loop  (1)  on  levels. 

(10)  Compute  Tfc,  k-1 ,K  hydrostatically  from  preliminary  estimates 

of  Z*  calculated  in  loop  (l)-(9).  Set  Q(T^)  equal  to  the  maximum 
of  Q(Z£)  «°d  Q(Z*+^).  If  either  of  Z£  or  Z*+1  is  missing,  then  T^  is 
missing. 

(11)  Calculate  final  Z£  for  k"l,K  hydrostatically  from  T^  and  Z^.  Now 

Z£  is  at  P£+j.  Zb  is  an  estimate  of  Z  at  the  bottom  of  layer  k.  For 
k“l,  Zb  •  Z^j  for  k>l,  Zb«  Z£_j  if  it  is  present  or  Z£_j  if  Zj^_j  is 
sdsslng.  If  T^  is  missing,  then  Z~  is  missing.  Quality  marks  are  set  so 
that  Q(Z* )  -  Q(Tk). 

(12)  Begin  loop  on  layers  k«i,  K.  This  loop  interpolates  (u,v)  and  In  q. 

Find  ft^,  the  first  t-l,  L  such  that  P^  <  P^.  Go  to  (15). 
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PTOSIG 


(13)  Fall  through  loop.  If  is  within  0.1  mb  of  Pk>  use  the  values  and 
quality  aarks  of  PL,  for  Ufct  Qfe  Q(Uk,  Vk)  and  Q(Qk>.  Go  to  (21). 
Otherwise, 

(14)  Do  not  interpolate.  Set  (uk>vk)>  Qk  to  missing.  Go  to  (21). 

(15)  If  £,-1  go  to  (14).  If  (P  ,  -  P  )  >  300  ab,  go  to  (14). 

*i  *i 

(16)  Find  lg,  the  first  level  at  or  below  £j-l  which  has  (U»V)  data.  Find  £u 

the  first  level  at  or  above  t^  which  has  (U,V)  data.  If  either  can't  be 

found  go  to  (19). 

(17)  If  (P.  -  P,  )  >  300  ab,  go  to  (19). 

*£  lu 

(18)  Interpolate  winds  at  and  £u  linearly  in  £n  p  to  Pk  to  obtain  (Uk»Vk). 

Set  Q(Uk,Vk)  to  maximum  of  (KU^.V^)  at  *•£  and  ®°  to  (20). 

(19)  Do  not  interpolate.  Set  (uk»vfc)  to  missing. 

(20)  Repeat  ( 16)— ( 19)  for  £n  (specific  hunidity).  (See  Mitchell,  1985.) 

(21)  End  of  loop  (12)  on  layers.  Return. 
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RE CALC 


Name:  RECALC  (subroutine) 

Purpose  :  Eliminates  the  observation  !  for  i  >  j  such  that  p^j,  the 

obaervation-observation  correlation  is  largest. 

Author:  D.  Norquist,  SASC,  1980  -  1986 


Documentation:  R.  Hoffman,  AER,  1986 


Referenced  by: 

ASAP2 

References: 

none 

Commons  used: 

none 

Arguments : 

RHO 

matrix  of  p^j 

RHS 

vector  of 

NNRHS 

pointer  to  pointer  arrrays 

NVRHS 

variable  types 

NEQS 

order  of  the  system 

AP 

work  space 

Description: 

(1)  Determine  position  of  maximum  p , j .  Only  the  row  position  imax  is 
actually  needed. 

(2)  Copy  arrays  to  temporaries ,  skipping  observation  i^^. 

(3)  Copy  temporaries  back  to  input  variables. 
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SATLTMP 


Name:  SATLTMP  (subroutine) 

Purpose:  Computes  sigma  layer  temperatures  and  sigma  level  heights 

from  satellite  observed  layer  temperatures,  using  the 
Flattery  algorithm  to  anchor  temperatures  and  first  guess 

**  a# 

or  Z£  to  anchor  heights*  (See  Norqulst,  1986b, 
pp.  35-36.) 


Author: 

D.  Norquis t ,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman,  AER, 

1986 

Referenced  by: 

MAST0R4 

v 

References : 

LOWTMPS 

Commons  used: 

DCONST,  FGDATA 

Arguments : 

P 

P£,  level  pressures 

T 

,  layer  temperatures 

ATA 

(A*A)“1  workspace 

NLVL 

L 

Zf 

TSIG 

PSIG 

0*0 

POH 

*0 

ZSTAR 

z* 

I/O  Units: 

none 

Description: 

(1)  Set  output  Z£  and  to  missing*  Return  if  any  of  the  input 
temperatures  (T  )  are  missing. 

X* 

(2)  Calculate  spaled  pressures,  s*  *  Pj/Pj  • 
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(3)  Use  Flattery  algorithm  (LOWTMPS)  to  obtain  level  temperatures 
(I*  ■  TSL)  from  the  layer  temperatures  (T  )• 

ft  ft 

(4)  Method  1  calculation  (not  presently  in  uae)  interpolates  to 

is  defined  as  missing  if  there  are  not  two  pressure  levels  surrounding 
lt»  or  if  one  of  the  surrounding  T£  values  is  missing.  (There  should  be 
no  missing  values  see  (1).)  Interpolation  in  linear  in  in  p.  Z~  is  then 
calculated  hydrostatically. 

(5)  Method  2,  which  used  400  A  as  the  prime  anchor  level,  starts  here.  Find 
£t,  the  first  1  such  that  p£  ■  p*  ,  trying  p*  -  400,  300,  200,  100,  SCO 
mb  in  succession.  If  none  found  return.  (Note  that  in  (1)  the  output  is 
preset  to  missing.) 

(6)  Integrate  downwards  from  level  ta  calculating  values  of  T£.  Of  the 
Flattery  temperatures  only  the  one  at  the  anchor  pressure.  It  is 
actually  used.  Given  the  upper  level  temperature  T*  and  the  layer  mean 

J L 

temperature  we  calculate  the  lower  level  temperature  for 

1  “  4a,  2,  -1.  The  calculation  assumes  T  is  linear  in  in  p  and 
that  T  is  a  mass  weighted  average.  (See  5.1.3.) 

ft 

(7)  Integrate  upwards.  Here  T£  and  x&  are  known  and  T£+i  is  calculated  for 
t  -  la,L-i. 

(8)  Interpolate  level  temperatures  T£  linear  in  in  p  to  obtain  sigma  layer 
values  T^.  xk  is  set  to  missing  if  TJ  values  above  and  below  it  are  not 
available  or  missing.  (There  should  be  no  missing  values;  see  (1)*) 

(9)  Determine  sigma  level  heights  Z£  by  intergrting  T^  hydrostatically. 

Original  anchor  for  height  integration  is  2*.  If  necessary  due  to  a 

0*0 

missing  T^  reanchor  to  first  guess  Z£.  In  practice,  missing  can  only 
occur  at  top  and  bottom  of  profile.  Normally,  then,  Z*  will  be  anchored 
to  Z*,  since  T^  will  usually  be  present. 
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Name:  SETFG  (subroutine) 

Purpose:  Initialises  the  blllnesr  Interpolation  by  reading  in  the 

southernmost  two  latitudes  of  the  first  quesa  grid*  Since 
data  are  sorted  by  latitude  only  two  latitudes  are  needed  at 
any  one  time* 


Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R«  Hoffman,  AER,  1986 


Referenced  by:  ASAP1 

References:  none 


Commons  used:  CCONST,  DCONST,  FGFLDS 


Arguments:  none 
I/O  units:  1*  5 
Description: 


(1)  Rewind  units  1,  5. 

M  M  M  M  ^ 

(2)  For  latitudes  1  and  2,  read  T,  U,  V,  R,  for  all  longitudes  and  layers 
from  unit  1.  Similarly,  read  Z*  for  all  longitudes  from  unit  5* 

Note:  File  structure  is  X,  k,  +,  while  /FGFLDS/  stores  data  by  X ,  +>  k. 
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4.  Input /Output  Usage* 

Program  ASAPI  uses  units  1-8,  10,  11,  INPUT  and  OUTPUT.  Table  4.1 
describes  these  files  briefly.  Table  4.2  lists  the  1/0  activity  in  ASAPI. 


Table  4.1 

Files  used  by  ASAPI.  All  files  are  either  formatted  (F)  or  unformatted  (U). 


Unit  Type 


Description 


1  U 


2  F 


3  U 


4  U 


5  U 


6  U 


7  U 


Fine  mesh  first  guess  values  for  interpolation.  Contains 
[<Xijk.  i«l,I,  k-i,K,  X-{T,U.V>), 
aijk.  i-l.I,  k-l,Kr),  (?nj,  i-l.I)].  J-l.J 
where  1*360,  J- 1 8 1 . 

Unpacked  GHE  Level  II  data.  See  Norqulst  (1984), 
pp.  17-33. 

Location  information  for  observations.  Contains: 

IN],  Un,  An,  daipJ  Bn,  points^],  n-l,N. 

Analysis  grid  first  guess  values.  Contains: 

[(Xijk,  k-l,K,  X*{T,U,V) ), 

<*ijk»  k-l,*r>*  **ijJ»  1-1 »*»  J-1*J 
where  1*61,  J*62. 

Fine  mesh  topography  for  interpolation.  Contains: 
lZ*tj,  i-l,I],  j-l,J 
where  1*360,  J-181. 

Analysis  grid  values.  Contains: 

((Xijk,  k-l,K,  X-{T,U,V>],  i-l,I,  j-l,J 
where  1*61,  J*62. 

The  sines  of  the  Gaussian  latitudes  (sin  4j],  j*l»J 
where  J*62. 
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8  U  The  residuals  used  by  the  analysis*  Contains: 

(OBD(n,NV),  tHI-1,73],  n-I,N. 

9  -  Not  used* 

10  U  If  this  Is  the  first  analysis  analysis  In  a  sequence,  this 

file  contains  the  NMC  BPS  at  the  Mandatory  levels  as  a 
function  of  analysis  grid  latitude: 

[<E  (XJt),  X  -  {Z,U,V})],  j-i,J 

where  >62,  L-12. 

Otherwise,  this  file  contains  the  BAB  for  height  fro*  the 
previous  analysis: 

k),  k"l,K]  ,  1*1,1, 
where  1*61,  J-62. 

11  U  The  EAE  and  Che  analysis  corrections: 

I(Ea(Xljk),  k-l.K,  X  -  (Z,U,V) ), 

<Xijk*  k"ltK»  *■<«*«»") >1»  J*1»J 

where  1-61,  J-62* 

INPUT  F  Standard  Input  device*  Contains  naeellst  DATINOP* 

OUTPUT  F  Standard  output  device.  Contains  printed  output. 

Note:  In  this  table  the  contents  of  each  record  Is  described  within  square 
brackets  (l  ]). 
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Table  4.2 


I/O  activity.  All  I/O  activity  is  listed  here  except  PRINTS  to  unit  OUTPUT. 


Balt  Activity  Occurences 

1  Rewind  SETFG  (1) 

Read  SETFG  (2) 

2  Open  AS  API  (3) 

iMd  ASAP1  (3),  (6),  (7),  (8) 

MASTOR1  (2),  (7) 

MASTOR2  (2),  (7).  (22) 

MAST0R4  (2),  (8),  (9),  (13),  (17) 
NAST0R6  (3) 


Rewind 

ASAP1  (7) 

Open-Close 

AS API  (3),  (6) 

MAST0R1  (18) 

MAST0R2  (23) 

MAST0R4  (18) 

Close 

MASTOR6  (19) 

3 

Write 

AS API  (11) 

4 

Read 

ASAP2  (1) 

5 

Rewind 

SETFG  (1) 

Read 

SETFG  (2) 

FC  (2) 

6 

Write 

ASAP2  (49) 

7 

Read 

ASAPl  (13) 

8 

Write 

AS API  (11) 

10 

Read 

ASAPl  (14) 

ASAP2  (4) 

11 

Write 

ASAP2  (49) 

Input 

Read 

ASAPl  (2) 

5.  Methods 

la  this  section  we  document  certain  methods  which  are  used  in  ASAP1. 

This  Section  of  documentation  Is  independent  of  the  actual  FORTRAN  code.  Some 
of  these  methods  are  used  repeatedly  in  the  code. 

5*1  Conversions 

5.1*1  Vector  wind  to  components 

Given  wind  speed  |v|  and  direction  0  the  eastward  and  northward 
wind  components  are  given  by 

u  -  -  |V|  sin  0 

V  »  -  j  V  |  C08  0 

Note  Chat  0  is  the  direction  the  wind  comes  from.  (See  also  Norquist  1986b, 
Hq.  23.) 

5.1.2  Heights  to  temperatures  hydrostatically  (and  vice  versa) 

The  integrated  form  of  the  hydrostatic  relationship  is 

Pu 

(g/R)  Az  -  £  T  din  p 

where  g  is  the  acceleration  of  gravity,  R  is  the  gas  constant  for  dry  air,  AZ 
is  the  thickness  of  the  layer  (Zu-Z4)  between  pressure  levels  p^  and  pu 
(Pg  >  pu),  T  is  the  temperature  and  p  is  pressure.  Strictly  T  is  virtual 
temperature.  If  T  is  linear  in  in  p,  then 

(g/R)  AZ  -  ¥  Ain  p 

where  T  is  layer  mean  temperature  is  given  by  (Tu  +  Tt)/2.  A  similar  equation 
is  given  by  Norquist,  (1986b,  Eq.  5). 

5.1.3  Layer  to  level  temperatures 

A  least  squares  procedure  for  converting  between  layer  and  level 
temperatures  due  to  Flattery  is  described  in  detail  by  Norquist  (1986b, 
Appendix  A,  pp  55-58.)  As  an  alternative,  if  one  level  temperature  is  known, 
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chen  all  other  level  temperatures  may  be  inferred  from  a  complete  set  of  layer 
temperatures  if  some  assumption  is  made  about  the  functional  form  of  T(p) 
between  the  levels.  If  T  is  linear  in  Jtn  p  then  one  obtains  the  method 
described  by  Norqulst  (1986b,  Appendix  D,  pp  64-65).  In  btief,  if  T  “  a  +  b 
In  p  between  p,  and  pu,  then  we  have 

pt  pt 

T  /  dp  -  /  (a  +  b  tn  p)  dp 

P  P 

*u  ru 

or 

(Pt-  PU)T  -  pt(a  +  b(in  p^-1))  -  Pu(«  +  b(tn  PU”D) 

Tu  -  a  +  b  tn  pu 

Tt  “  a  +  b  in  P( 

which  are  three  equations  in  the  three  unknowns  a,  b  and  either  Tu  or 
E.g.,  if  Tu  is  known,  then 

a  -  Tu  -  b  in  pu 

b  -  (T  -  Tu)  (Pjl-  PU)/(PU  ♦  P1Un(pjt/pu)-lJ) 

5.1.4  Humidity 

Dew  point,  T^,  dew  point  depression  T-T^,  specific  humidity  q  and 
relative  humidity  r  are  all  used  in  the  analysis.  These  are  related  by  the 
following  expressions.  (Dutton,  1976,  Chap.  8): 

r  -  e/es 

e8  "  eo  eXP  ^ R  T ^  ^ 

V  O 

-  6.11  exp(19. 9274-5443. 3618/TJ 
e  -  ea  (Td) 

W  -  (1/e)  e/(p-e);  e  -  eWp/(i+eW) 
q  -  W/(l+W);  W  -  q/(l-q) 
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where 

e  “  (molecular  weight  air )/(oolecular  weight  water)  *  1.61, 

e  is  the  vapor  preaaure,  eg  the  aaturation  pressure,  and  W  the  mixing  ratio 
(i.e.  the  ratio  of  the  mass  of  vapor  to  the  mass  of  dry  air).  Similar 
equations  are  given  by  Norquist  (1986b,  Eq.  1-4  and  22)  except  that  volume 
mixing  ratio  is  used  and  the  definition  e  is  reversed  in  Eq.  1-4. 

5.2  Generic  quality  marks 

Wd  dletinqulsh  (a)  quality  indicators  Qj,  which  are  read  from  the  GWE 
Level  11  data,  (b)  generic  quality  marks,  Q,  which  are  used  in  the  analysis 
selection  procedures  and  (c)  quality  levels,  QL,  which  are  used  in  the  buddy 
check  procedures. 

The  quality  indicators  Qj  are  described  by  WHO  (1978),  e.g.  Table  IV  of 
Appendix  A.  Ordinarily,  Qj  “  0  indicates  no  quality  control  check  (QCC)  was 
made  and  Qj  -  1,2...  indicate  decreasing  quality  was  determined  during  the 
QCC.  However  Qj  has  different  meanings  for  different  data.  In  ASAP1  all  Qj 
are  translated  to  Q,  which  always  have  the  same  meaning  (see  Table  5.1).  Note 
If  Q(Z1)  <  Q(Zj )  then  Is  better  than  Zy 
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Table  5.1  Generic  Quality  Mark* 

Quality  Mark  _ Meaning 

0  Value  found  correct  during  QCC. 

1  No  QCC  was  made* 

2  Value  found  erroneous  during  QCC,  reconstituted  value 
Inserted* 

3  Value  found  missing  during  QCC,  reconstituted  value  inserted* 

4  Value  found  missing  during  QCC,  new  value  assigned* 

5  QCC  made,  new  value  probably  entered* 

6  Value  found  suspect  during  QCC* 

7  Value  found  erroneous  during  limits  checks* 

8  Value  found  erroneous  during  QCC* 

9  Value  missing. 


The  quality  levels  QL  used  In  the  buddy  check  are  given  by  Norquist 
(1986b,  pp  7—8). 

1Q+1 

2(0+3)  +  1 
3(Q+3)  +  2 
3(Q+3)  +  3 

5.3  Buddy  Check  Boxes 

The  following  procedure  determines  the  buddy  check  box  index  Bn,  given 
the  latitude  6  and  longitude  X  of  observation  (Norquist,  1986b,  p.5): 

(1)  Determine  the  latitude  band  i  *  1,...  18  from  SP  +  NP  in  10*  increments 
by  i  -  ($  +  90)/ 10  +  1.  The  center  of  the  latitude  band  is  then 
♦c  -  101  -  95. 

(2)  Determine  the  number  of  boxes  for  this  band  as  n^  ■  36  *  cos  4>c . 
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(3)  Determine  the  longtitude  Index. 

j  ■  A/AX  +  1  where  AX  “  360/n£  . 

(4)  Set  Bn  -  1001  +  J. 

Mote  that*  for  the  CRAY-1,  X  and  4  oust  be  Incremented  by  0*001  to  get 
reproducible  results.  This  should  have  no  effect  since  X  and  $  are  coded  to 
hundredths  only* 


5.4  Vertical  grid  and  interpolation 

A  sigma  coordinate  system  is  used.  The  interface  value  o*  are  given  by 
p*/p^  for  k“l,  K+l  .  The  layer  thicknesses  Aofc  ■  o£  -  o£+1  are  defined  as 
0*075,  0*125,  0*150,  0.150,  0. 125,'' 0.075,  6*0.050.  In  addition  an  assumed 
underground  layer  has  Aou  “  0*075*  Indexing  starts  with  the  ground  k-1  and 
ends  with  the  top  of  the  atmosphere  k  -  K+l.  The  layer  ok  values  are  given  by 

o^^J/CCx+DAo^)]1^ 

where  x  ■  R/Cp  (Brenner  et  al*,  1982,  Eq.  32;  Norqulst,  1986b,  Fig.  1). 

Note  that  the  vertical  indexing  for  height  is  generally  offset  by  1. 

Thus,  Z£  is  the  height  at  • 

Vertical  interpolation  is  generally  linear  in  inp;  with  data  X£  and  3^  at 
PA  and  pu  we  have 

X-X£  tn(p/p£) 

Vxt "  *■<*>,> 

In  ASAP1,  this  relationship  is  often  put  in  point-slope  form  using  the 
midpoint  of  the  interval, 

X-(X£+  Xu)/2  +  [<X£-  Xu)/tn(p£/pu)]  (In  p  -  1/2  In  p£pu) 
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5.5  Polar  stereographic  projection 

For  gridpoints  with  |<j>  |  >  70*  a  polar  stereographic  projection  is  used 

8 

for  data  selection,  buddy  checking  and  performing  the  analysis*  Nap  factors 
are  used  throughout  following  Bergman  (1979)  although  Dey  and  Morone  (1985) 
suggest  that  for  a  polar  cap  the  map  factors  are  not  necessary* 

The  conversion  of  (X,$)  coordinates  to  (x,y)  coordinates,  in  a  system  in 
which  the  positive  x  direction  points  from  the  pole  towards  Greenwich,  is 
given  by  (Bergman,  1979,  Eqs.  3. 3-3.4) 

m  »  2/(1  +  sin  (±  $)). 


x  ■  am  cos  $  cos  X 


y  -  *  am  cos  $  sin  X 

where  m  is  the  map  factor,  a  is  the  earth's  radius  and  the  +  (or  -) 

corresponds  to  the  Northern  (or  Southern)  Hemisphere.  The  distance  djj 
between  two  points  indexed  i  and  j  is  jthen  given  by  (Bergman,  1979,  Eq.  3.5) 


i-  [(xi-  Xj)2  +  (y£-  y^  )2 J 


where 


m 


(mi  +  “j)/2* 


In  the  coordinate  system  chosen  for  the  projection,  only  the  winds  at 
90°W  are  unchanged  from  the  (X,$)  system.  The  others  must  be  rotated  by 
(Norqulst,  1986b,  Eq.  6) 


/  sinX 
\*cosX 


±C08X 

sinX 


where  (u',v*)  is  the  wind  in  the  stereographic  projection  and  (u,v)  is  the 
wind  in  the  (X,$)  system.  Again  +  (or  -)  corresponds  to  the  Northern  (or 
Southern)  hemisphere.  The  reverse  conversion  is  given  by  the  transpose  of  Che 
above  rotation  matrix: 
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/  ainA 

\  ±C08X 


♦coeA 

sinA 


) 


5*6  Coefficient  of  geoetrophy 

•  The  coefficient  of  geoetrophy  G  decouples  the  wind  and  height  analysis  in 

the  tropica*  Bergman  (1979,  Eq.  3.20)  uses 


G  -  1  -  exp  (U|/20) 


while  Dey  and  Norone  (1985,  Eq.  BIO)  use 


G 


0  0<|*|<10 

[cos  (12 1 4  | -300)  +  l)/2  10<|+|<25 

1  25<  | <  90 


AS API  uses  (Norquist,  1986b,  pl4,  and  Fig.  2) 

G  (♦)  G  (25)  Ot  |a  |<25 

G  m 

G  (4)  25<U|<90 
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